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CHAPTER  I  INTRODUCTION 


1.1  STIFF  EQUATIONS 

Systems  of  ordinary  differential  equations  which  arise  in  the 
description  of  physical  problems  often  have  greatly  differing  time  con¬ 
stants  and  are  very  stable  (Bjurel  et  al.,  1970).  They  are  known  as  stiff 
equations.  A  model  stiff  equation  is  given  by 

y' (t)  =  Ay(t)  y(tQ)  =  n 

where  yelRm,  A  is  an  mxm  real  matrix  with  all  its  eigenvalues  A-|,A2,. •  •  »xm 

lying  in  the  open  left  half  complex  plane,  and  the  ratio 

| Re  A  .  | 

max  - 

i.j  | Re  A. | 

is  large.  All  solutions  y(t)  =  exp[(t-tQ)A]n  approach  zero  as  t  increases, 
hence  the  system  is  asymptotically  stable  about  y  =  0.  The  component  of 
y(t)  corresponding  to  an  eigenvalue  a.  whose  real  part  is  large  in  magnitu- 

J 

de  will  soon  become  negligible  relative  to  the  dominant  component,  i.e., 
the  one  corresponding  to  the  eigenvalue  whose  real  part  is  the  smallest  in 
magnitude,  which  could  still  be  large  enough  to  be  of  interest.  For  exam¬ 
ple,  we  consider  the  following  equation  when  m  =  2: 


A-l 

-A-l 

y^t) 

y,U0) 

y'2(t) 

-A-l 

A-l  1 

J 

y2(t) 

*2  (*())_ 

_n2_ 

where  A  «  0.  The  true  solution  is 


-t  At 

y-|  (t)  =  c-| e  +  c2e 

O 

— j 

II 

zi 

+ 

Z3 

ro 

ro 

-t  At 

y2(t)  =  c-j e  -  c2e 

c2  =  (n-|  -  n2)/2 
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The  component  exp(xt)  is  negligible  after  t  >_  1 ,  and  thus  we  would  like  to 
be  able  to  use  a  larger  stepsize  h  from  that  time  on  when  integrating  the 
system  numerically,  since  we  need  only  follow  the  dominant  component 
exp(-t).  However,  numerical  stability  of  methods  which  are  not  specifical¬ 
ly  designed  for  stiff  equations,  such  as  the  one-step  Euler's  method,  re¬ 
quires  that  the  value,  for  the  m-dimensional  case, 

max  |hxJ 
l<j<m 

remain  small  throughout  the  range  of  integration,  forcing  the  stepsize  h  to 
be  unnecessarily  small  (Lambert,  1973,  pp. 229-231). 

It  would  then  obviously  be  desirable  to  have  methods  which  do  not 
require  |  hx .  |  to  be  small  except  for  accuracy.  Such  relaxation  of  the  sta- 
bility  requirement  is  the  essential  property  of  methods  for  stiff  problems. 

In  general,  we  are  concerned  with  solving  initial  value  problems 
involving  systems  of  first  order  ordinary  differential  equations  (ODE)  of 
the  form 

y'(t)  =  f(t,y(t))  y(t0)  =  n  te[t0,T]  (1.1) 

where  f  is  continuous  in  t  and  satisfies  a  Lipschitz  condition 

I  f(t,y-j )  -  f ( t ,y2 ) |  i  L|y1  -  y2l 

for  some  positive  constant  L  in  the  region  tQ  ±  t  <1,  |y|  <  °°,  so  that  a 
unique  solution  y(t)  exists  and  is  in  C^[tg,T]  (Coddington  and  Levinson, 
1955,  p.10).  We  further  assume  that  f  has  continuous  partial  derivatives 
with  respect  to  y  and  that  the  solution  y  has  as  many  continuous  deriva¬ 
tives  as  is  necessary. 
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1.2  PROPERTIES  OF  NUMERICAL  METHODS  FOR  STIFF  SYSTEMS 


Most  numerical  methods  for  integrating  the  initial  value  problem 
(1.1)  approximate  the  true  solution  y(t)  at  a  sequence  of  mesh  points  tQ  < 
t-j  <  ...  <  tN  =  T.  We  denote  the  numerical  approximation  to  y(tn)  by  yn,  0 
<_  n  £  N.  The  difference  between  successive  mesh  points,  tn  -  tn_-j  •  is  cal¬ 
led  a  stepsize,  which  is  not  necessarily  the  same  for  different  steps. 
However,  for  ease  in  analysis,  we  assume  that  the  stepsize  is  constant 
throughout  the  range  of  integration. 

When  using  a  numerical  method,  we  are  concerned  with  how  accurate 
we  can  make  the  numerical  approximations  to  the  true  solution  by  picking 
one  or  more  parameters,  for  example,  the  stepsize.  A  method  is  said  to  be 
convergent  if  for  any  initial  value  problem  (1.1), 


lim  max  |y  -  y(t  )|  =  0 
N-*=°  0<n<N 


in  which  t^  =  T.  Thus,  when  using  a  convergent  method  to  integrate  (1.1), 
we  can  achieve  any  desired  accuracy  by  picking  a  small  enough  stepsize  h. 
Another  concept  concerning  error  propagation  is  also  important  since,  in 
practice,  computations  are  hardly  exact  because  of  the  finite  number  of  di¬ 
gits  that  can  be  carried.  We  loosely  define  a  method  to  be  stable  if,  for 
each  problem  of  the  form  (1.1),  there  exists  a  hQ  >  0  such  that  a  perturba¬ 
tion  in  the  equations  defining  yn  by  a  fixed  amount  produces  a  bounded 
change  in  the  numerical  solution  when  any  he(0,hg]  is  used  as  stepsize.  A 
more  precise  definition  of  stability  is  given  in  Stetter  (1973,  p.9).  Both 
the  concepts  of  convergence  and  stability,  necessary  for  any  method  to  be 
useful,  are  concerned  with  the  limiting  process  as  h  +  0. 

In  practice,  we  are  more  concerned  with  the  effect  of  error  accu¬ 
mulation  --  not  only  roundoff  error  as  introduced  earlier,  but  also  discre- 
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tization  error  from  approximating  the  true  solution  --  when  some  finite 
stepsize  h  is  used.  We  thus  introduce  the  concept  of  absolute  stability 
for  a  given  fixed  stepsize.  Since  the  concept  is  problem  dependent,  we  de¬ 
fine  absolute  stability  for  the  class  of  differential  equation 

y'(t)  =  xy(t)  y(tQ)  =  n  (1.2) 

where  x  is  a  complex  constant. 

Definition 

The  region  of  absolute  stability  4  of  a  method  is  the  set  of  all 
hx  such  that  when  the  method  is  applied  with  stepsize  h  >  0  to  equation 
(1.2),  the  numerical  solution  yn  ->  0  as  n  -»•  ®. 

The  method  is  then  absolutely  stable  for  stepsize  h  and  for  the 
equation  (1.2)  if  hxeA.  The  following  terminology  has  been  used  to  descri¬ 
be  methods  according  to  the  extent  of  their  region  of  absolute  stability: 

Definition  (Cryer,  1973) 

A  method  is  AQ-stable  if  the  negative  real  axis  (z  :  Im  z  =  0, 

Re  z  <  0}  is  in  A. 

Definition  (Widlund,  1967) 

A  method  is  A(a)-stable,  ae(0,V2),  if  the  wedge  Sa  =  (z  : 

|Arg ( -z ) |  <  a,  z  f  0}  is  in  A,  The  a  angle  of  the  largest  such  wedge  is 
called  the  angle  of  absolute  stability  associated  with  the  method.  A 
method  is  A(Q)-stable  if  it  is  A(os)-stable  for  some  ae(0,ir/2),  and  is 
A(tt/2) -stable  if  it  is  A(a)-stable  for  all  ae(0,V2). 
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Definition  (Dahlquist,  1963) 

A  method  is  A-stable  if  its  region  of  absolute  stability  contains 
the  open  left  half  complex  plane. 

From  the  definition,  it  is  obvious  that  A-stability  is  equivalent 
to  A(Tr/2)-stability.  Moreover, 

A-stability  =>  A(a)-stability  A(0)-stabil ity  =>  A0-stability 

Methods  which  are  Ag-stable  shall  be  loosely  called  stiff  methods. 

The  angle  of  absolute  stability  is  only  one  of  a  number  of 
parameters  which  have  been  proposed  for  measuring  the  extent  of  the  region 
of  absolute  stability  A  (cf.  Gear,  1971,  p. 21 3 ,  Odeh  and  Liniger,  1971, 
Gupta,  1976).  But  it  is  probably  the  best  such  measure,  especially  for 
methods  with  automatic  stepsize  selection.  When  such  methods  are  applied 
to  the  equation 

y'(t)  =  xy (t)  +  g(t)  x  complex 

the  integration  starts  out  with  hx  near  the  origin,  and  as  the  integration 
proceeds,  the  stepsize  h  increases  and  the  value  hx  moves  away  from  the 
origin.  However  if  hx  approaches  the  boundary  of  A,  this  would  be  detected 
by  the  error  estimator,  and  any  further  movement  of  hx  would  be  prevented. 
The  implication  of  this  is  that  not  all  of  A  is  "used",  but  only  that  por¬ 
tion  which  can  be  reached  from  the  origin  by  rays  lying  entirely  inside  A, 
The  most  "desirable"  portion,  where  any  h  >  0  can  be  used,  is  described  by 
the  angle  of  absolute  stability.  This  parameter,  a,  also  serves  as  a  good 
indicator  of  the  problem  class  for  which  the  method  is  suitable,  namely, 
those  problems  for  which  the  large  eigenvalues  of  the  Jacobian  lie  inside 
the  wedge  S  .  Hence,  in  an  ODE  solver,  a  can  be  an  extra  parameter  which 
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is  used  to  identify  among  a  family  of  methods  of  order  k  the  A(a) -stable 
method  that  should  be  used.  The  value  of  a  can  be  either  supplied  by  the 
user  or  chosen  automatically  by  some  detecting  device  in  the  code.  It 
would  then  be  desirable  to  have  A(a)-stable  methods  of  any  order  for  any 
o«s(0,ir/2]. 

Since  the  true  solution  of  (1.2)  for  any  x  lying  in  the  open  left 
half  plane  converges  to  zero  as  t  +  °°,we  would  prefer  methods  whose  region 
of  absolute  stability  contains  the  open  left  half  plane,  i.e.,  A-stable,  so 
that  the  numerical  solution  yn  also  converges  to  zero  as  n  ~  for  any 
stepsize  h  >  0.  On  the  other  hand,  if  X  is  not  in  the  open  left  half 
plane,  y(t)  will  not  converge  to  zero,  and  in  fact  will  diverge  if  Re  X  > 

0.  However,  if  hXGA,  yn  converges  to  zero  and  consequently  invalid  solu¬ 
tions  are  produced  (Lindberg,  1974).  Hence,  methods  whose  region  of  abso¬ 
lute  stability  is  equal  to  the  open  left  half  plane  are  desirable. 

We  end  this  section  with  another  definition: 

Definition  (Odeh  and  Liniger,  1971) 

A  method  is  AOT-stable  if  its  region  of  absolute  stability  con¬ 
tains  a  neighborhood  of  infinity. 

A^-stability  is  desirable  especially  for  those  problems  which 
have  a  slowly  varying  component  and  a  rapidly  decaying  component,  so  that 
transient  that  decays  to  zero  rapidly  in  the  true  solution  would  not  be 
present  in  the  numerical  solution  as  slowly  dampened  oscillatory  com¬ 
ponents. 

1.3  LINEAR  MULTISTEP  METHODS 

We  restrict  our  study  to  the  class  of  linear  k-step  formulas 
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k 

E 

0=0 


ajyn+j 


k 

h  E 
0=0 


pjf(WW  = 


n  =  0,1 .... 


(1.3) 


where  k  is  some  fixed  positive  integer,  a,,  e.,  j  =  0,1,.. .,k,  are  given 

J  J 

constants  such  that 


ak  * 0 

l“0l  +  I eo I  ^  ^ 


and  h  >  0  is  the  stepsize  assumed  to  be  fixed  throughout  the  range  of  in¬ 
tegration. 

When  k  >  1,  the  first  step  of  integration  requires  numerical 

values  for  yQ,y1 . yk1-  Since  we  are  only  given  n  as  an  initial  value, 

the  k  values  have  to  be  computed  from  i>  A  common  technique  is  to  use  a 
one-step  method  such  as  the  Runge-Kutta  method  (Gear,  1971,  Chapter  2)  to 

compute  yQ,y1 . yk_r  The  procedure  used  to  compute  yQ,y1 ,. . . ,y|<_1  from 

n  is  called  a  starting  procedure.  The  k-step  formula  (1.3)  together  with  a 
starting  procedure  constitutes  a  k-step  method.  The  starting  procedure  is 
said  to  be  consistent  with  (1.1)  if  for  j  =  0,1,..., k-1, 

yj  =  yj(n,h)  n 


as  h  ->  0. 

At  each  step  of  the  integration,  in  the  case  when  ek  f  0,  i.e., 
when  the  formula  is  implicit,  we  have  to  solve  a  system  of  nonlinear  equa¬ 
tions  for  yp+k.  For  stiff  equations,  this  is  usually  accomplished  by  some 
kind  of  modified  Newton  iteration.  Results  are  usually  satisfactory  if  the 
stepsize  h  is  small  enough  and  if  the  initial  estimate  is  sufficiently  ac¬ 
curate  (Lambert,  1973,  p.239).  The  initial  estimate  could  be  obtained  by 
using  an  explicit  formula. 
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Unless  explicitly  specified,  the  following  assumptions  about  the 
linear  k-step  formula  (1.3)  are  made: 
k 

(A)  z  3.  =  1,  a  normalization.  (1.4) 

j=0  J 

(B)  The  formula  (1.3)  satisfies  the  root  condition,  i.e.,  all  the 
roots  of  the  polynomial 

k  j 

p(0  =  z  <*.e  (1.5) 

j=0  3 

lie  inside  the  unit  circle  |c|  1  1  and  those  on  |c|  =  1  are  sim¬ 
ple.  Occasionally,  we  require  that  (1.3)  satisfies  the  strict 
root  condition,  i.e.,  all  roots  of  p(c)  lie  inside  |c|  <  1  except 
for  a  root  at  s  =  1 . 

(C)  The  order  of  the  formula  (1.3)  is  k,  where  formula  (1.3)  is  said 
to  be  of  order  p  if  cQ  =  c^  =  ....*  c  s  0  in  the  following  Taylor 
series 

k 

e  [o1y(t+jh)  -  he.y'(t+jh)]  = 

0=0  J  3 

=  cQy(t)  +  c-j  hy 1  (t)  +  ...  +  cphpy^(t)  + 

+  cp+]hP+VP+1)(t)  +  ...  (1.6) 

It  is  said  to  be  of  exact  order  p  if  cQ  =  c^  =  ...  =  cp  =  0  and 
cp+-j  4  0,  in  which  case,  we  call  cp+^  the  error  constant  of  the 
formula.  Formulas  of  order  at  least  1  are  said  to  be  consistent. 

(D)  The  polynomials  p(c)  and  cr(c)  have  no  common  factors,  where  p(c) 
is  as  defined  in  (1.5)  and 
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k  i 

c(c)  =  z  e,?J  (1.7) 

j=o  J 

Stetter  (1973,  pp. 186-1 87)  shows  that  if  the  polynomials  p  and  a 
of  a  linear  multi  step  formula  possess  one  or  more  common  factors, 
there  exists  an  equivalent  formula  with  a  smaller  step  number  in 
the  sense  that  they  generate  the  same  numerical  solution  from 
given  starting  values. 

We  also  assume  that  the  starting  procedure  is  consistent.  The  k-step 
method  is  said  to  be  consistent  if  both  the  starting  procedure  and  the  k- 
step  formula  are  consistent.  It  is  known  that  a  k-step  method  is  stable  if 
and  only  if  the  k-step  formula  satisfies  the  root  condition  (Stetter,  1973, 
p.208).  It  has  also  been  proved  that  a  k-step  method  is  convergent  if  and 
only  if  it  is  stable  and  consistent  (Stetter,  1973,  pp.13,  211-213  --  note 
that  some  of  the  details  are  missing). 

Many  authors  require  that  the  formula  (1.3)  together  with  an  ar¬ 
bitrary  consistent  starting  procedure  be  convergent  (e.g.,  Henrici,  1962, 
p.218).  However,  Stetter  (1973,  p.78)  indicates  that,  except  in  pathologi¬ 
cal  cases,  this  requirement  is  no  stronger  than  ours. 

It  should  also  be  remarked  that  the  assumptions  about  constant 
stepsize  and  order  are  just  idealizations  that  make  our  analysis  tractable. 
In  practice,  both  the  stepsize  and  order  will  vary.  Also,  in  practice, 
formulas  that  satisfy  the  strict  root  condition  are  more  desirable.  We 
call  those  formulas  strongly  stable. 

We  denote  a  k-step  method  by  (p,a)  where  p  is  as  defined  in  (1.5) 
and  a  in  (1.7),  and  use  the  symbol  L[k]  to  represent  the  class  of  (p,ct) 
which  satisfy  conditions  (A)-(D).  It  is  easily  seen  that,  under  assump- 
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tions  (A)  and  (C),  we  have 


c(l)  =  l 

cQ  =  p(D  =  0 

c1  =  p'(l)  -  0(1)  =  0 


and,  in  general,  for  q  =  2,3,... 


1 


s  (j'V-  -  qjq_13,) 


H  q!  j=l 

VJe  remark  that  in  order  for  a  method  (p,a)  in  £[k]  to  be  Aq- 
stable,  it  is  necessary  that  (Cryer,  1973) 


6kt  0 

order  £  k  except  the  trapezoidal  rule 
Furthermore,  the  region  of  absolute  stability  A  associated  with  a  method 
(p,a)  can  be  described  by 

A  =  (yeC  :  |tj(y)|  <  1,  j  =  l,2,...,k> 


where  c j (p ) ,  j  =  1,2 . k,  are  the  k  roots  of  p(t)-ycr(s),  depending  on  the 

value  of  yeC  (Lambert,  1973,  p.222). 

The  stability  condition  requires  that  all  the  roots  of  a  po¬ 
lynomial  lie  inside  the  unit  circle.  Since  there  exists  a  reasonably  sim¬ 
ple  algebraic  condition  for  all  the  roots  to  lie  in  the  left  half  complex 
plane,  we  wish  to  transform  the  unit  circle  to  the  left  half  plane.  This 
can  be  done  by  using  the  bilinear  transformation 

z+1 

£=—  (1.8) 

z-1 

which  maps  the  left  half  plane  in  complex  z-space  1-1  onto  the  unit  circle 
in  complex  s-space.  Apply  the  transformation  (1.8)  to  the  polynomials  p(?) 
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and  a(c)»  and  define  the  corresponding  polynomials  in  z  by 

k  .  .  z+1 

r(z)  =  2  a.zJ  =  (z-l)kp( - ) 

j=0  3  z-1 

k  .  .  z+1 

s  (z )  =  z  b.zJ  =  (z-l)ka( - ) 

j=0  3  z-1 

It  is  easily  seen  that 

ak  =  p(l)  =  0 

bk  =  <t(1)  =  1 

Furthermore,  a  method  is  of  order  k  if,  and  only  if,  the  following  holds: 


z+1  2k+  c, 


- log  - 

s(z)  z-1 


_ 5li  +  o( _ ) 

zk+l  zk+2' 


as  z  ^  oo 


(obtained  from  the  definition  of  order  by  setting  y(t)  =  exp(t)  in  (1.6)). 

Therefore,  requiring  that  the  method  be  of  order  k  imposes  k  additional 

linear  conditions  on  the  a.'s  and  the  b.'s: 

3  J 


k  b  • 

a-  =  2  I  -  j  =  0,1 ,... ,  k-1  (1.9) 

J  i=j+l  1-j 
i-j  odd 

Hence,  the  method  is  uniquely  parameterized  by  the  k  parameters 
b0,b1,...,b(<_-|  (bk  =  1).  Occasionally,  we  denote  such  parameterization  by 

(bQ,b-| . bk-i)»  or  by  tbe  polynomial  s(z),  or  by  (r,s)  in  analogy  to 

(p.o). 

The  parameterizations  (p,a)  and  (r,s)  are  related  by  the  fol¬ 
lowing  linear  transformation  (Duff in,  1969): 

a  =  Met  b  =  MB 


where  a  =  [aQ,a-|,...,akj  ,  b  =  [bQ,b-|,...  »bk3  *  a  -  [ao*al*’  *  *  ,ak^  *  ^ 
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[3q , B-j  , 3j<]T ,  and  M  =  [M. .]  is  a  (k+l)x(k+l)  matrix  whose  elements  are 

■  J 

given  by 


M10  =  (-1 )k_1 (  ) 

Mkj  =  1 

Mij  "  +  Mi+l,J-l  +  K1+l,j 


with  M’1  =  2'kM. 


i  =  0,1,..., k 
j  -  0,1 . k 

i  =  0,1 . k-1  j  =  1 ,2,. . .  ,k 


We  note  that  the  polynomial  r(z)  of  a  stable  method  has  roots  all 
in  the  left  half  complex  plane,  with  at  most  simple  roots  on  the  imaginary 
axis,  and  likewise  for  the  polynomial  s(z)  of  an  A(0)-stable  method 
(Widlund,  1967)  --  the  condition  is  relaxed  for  Ag-stable  methods  whose 
s(z)  could  have  at  most  double  roots  on  the  imaginary  axis  (Cryer,  1973). 
The  case  when  all  the  roots  of  r(z)  (  s(z)  )  lie  in  the  open  left  half 
plane  corresponds  to  the  method  being  strongly  stable  (  A^-stable  ). 

The  parameterization  (r,s)  has  several  advantages  over  (p,a). 

The  restriction  that  (r,s)  be  of  order  k  imposes  a  simple  relation  on  the 
a.'s  and  the  b^.'s  (cf.  (1.9)).  Moreover,  the  error  constant  c^  can  be 
expressed  simply  as  a  linear  combination  of  the  even  b-'s  (cf.  (2.3)). 

J 

Also,  the  Hurwitz  criterion  (cf.  Appendix  A)  for  polynomials  with  roots  in 
the  open  left  half  plane  is  "simpler"  than  the  Schur  criterion  (Marden, 
1966,  p. 1 98)  for  polynomials  having  roots  inside  the  unit  circle.  However, 
the  aj's  and  the  bj's  have  no  natural  interpretation  like  the  cu's  and  the 
3j's  have,  namely,  parameters  of  the  linear  difference  equation  defining 
the  numerical  solution  yn. 

We  finally  remark  that  several  versions  of  the  parameterization 
(r,s)  have  been  used  by  different  authors.  The  one  described  in  this  sec- 
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tion  is  used  in  Liniger  (1975),  Genin  (1973),  and  Jeltsch  (1976). 

Dahlquist  (1963),  Widlund  (1967),  and  Cryer  (1973)  use  the  same  bilinear 

transformation  (1.8)  but  define  the  polynomials  r(z)  and  s(z)  differently: 

z-1  L  z+1 
r(z)  =  (— )  p(— ) 

2  z-1 

z-1  k  z+1 
s(z)  =  (—  )o(— ) 

2  z-1 


Henrici  (1962,  p. 230 )  and  Gear  (1971,  p. 1 95 )  use  a  different  bilinear 

transformation  and  hence  different  r(z),  s(z): 

1+z 
s  =  — 

1-Z 


k 


r(z)  =  ( — )np( 
2 


1-z  k 

s(z)  -  ( — )  °( 

2 


1+z 

— ) 
1-z 

1+z 
— ) 
1-z 


1.4  OBJECTIVE 

As  pointed  out  in  section  1.2,  in  order  to  be  able  to  solve  gen¬ 
eral  stiff  systems  of  the  form  (1.1),  we  would  very  much  prefer  to  use 
methods  which  are  A-stable,  so  that  they  would  be  suitable  for  all  stiff 
problems.  However,  Dahlquist  (1963)  proves  that  the  maximum  order  for  A- 
stable  linear  multi  step  methods  is  only  2.  To  allow  for  methods  having  or¬ 
der  greater  than  2,  we  have  to  relax  our  stability  requirement  and  consider 
A(Q)-stable  methods.  Cryer  (1973)  shows  that  there  exist  Ag-stable  methods 
of  arbitrary  order  k  by  explicitly  constructing  a  class  of  Ag-stable 
methods  in  £[k],  which  are  later  shown  by  Jeltsch  (1976)  to  be  in  fact 
A(0)-stable.  However,  the  error  constant,  which  is  an  asymptotic  measure 
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of  the  accuracy  of  the  method  (cf.  Chapter  II),  of  Cryer's  methods  becomes 
astronomical  in  magnitude  when  k  is  large,  say  when  k  =  10.  Consequently, 
we  are  compelled  to  find  more  accurate  A(0)-stable  methods. 

Work  has  been  done  in  seeking  stiff  methods  with  as  high  order  as 
possible.  Widlund  (1967)  shows  that  A(a)-stable  methods  exist  for  any 
a£(0,Tr/2)  for  the  cases  when  k  =  3,4.  Using  an  interactive  computer  gra¬ 
phics  program.  Dill  and  Gear  (1971)  find  stiffly  stable  methods  (definition 
in  Gear  (1971,  p . 21 3 )  --  note  that  stiff  stability  implies  A(0)-stabil ity) 
of  orders  7  and  8.  They  only  consider  the  class  of  methods  most  of  whose 
p.'s  being  zero.  However,  Skeel  (1977)  shows  that  every  k-step  method  is  a 

J 

(k+l)-value  method.  The  same  observation  is  mentioned  in  Wallace  and  Gupta 

(1973)  and  is  proved,  for  the  special  case  when  the  method  is  of  order  k+1 , 

in  Osborne  (1966).  Hence  insisting  that  most  of  the  p.'s  be  zero  would  not 

J 

reduce  the  number  of  previous  values  needed  to  be  stored.  Jain  and  Srivas- 
tava  (1970)  consider  the  polynomial 

a(c)  =  ek"r(?-c)r  0<r<k,  ce[-l,l) 
and  find  stiffly  stable  methods  of  order  as  high  as  11  for  appropriate 
choices  of  r  and  c.  Gupta  and  Wallace  (1975)  parameterize  a  k-step  method 
by  what  they  call  a  modifier  polynomial  C(x).  By  requiring  that  C(x)  or 
C'(x)  approximates  zero  in  some  sense  (e.g.  interpolation,  Chebyshev,  least 
squares),  they  succeed  in  finding  stiff  methods  of  order  up  to  9  with  a  = 
71.4°  (the  FMPD60  formulas  in  Gupta  (1976)). 

The  aim  of  our  study  is  to  discover  some  of  the  limitations  on 
the  accuracy  of  stiff  methods  in  L[k],  We  first  choose  a  measure  of  accu¬ 
racy  for  our  study,  then  we  investigate  the  following  problems  concerning 
the  class  of  stiff  methods  in  L[k]: 
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(i)  How  accurate  can  they  be  ? 

(ii)  For  a  given  ae(0,V2),  what  is  the  limitation  on  the  accuracy  of 

A(a)-stable  methods  ? 

(iii)  For  a  given  constant  C  >  0,  what  is  the  limitation  on  the  angle 
of  absolute  stability  for  methods  having  error  constant  of  mag¬ 
nitude  C  ? 

(iv)  For  any  cus(0,ir/2),  do  there  exist  methods  which  are  A(a)-stable  ? 

If  they  exist,  how  accurate  are  they  ? 

Chapter  II  is  devoted  to  a  discussion  on  measures  of  accuracy. 

The  answers  for  the  above  questions  for  the  cases  when  k  =  1 ,2  are  given  in 
Chapter  III,  where  questions  (i)  and  (iv)  for  a  general  k  are  also  con¬ 
sidered.  An  algorithm  for  solving  minimax  problems  numerically  is  introdu¬ 
ced  in  Chapter  IV.  Questions  (ii)  and  (iii)  for  a  general  k  are  formulated 
as  a  minimax  problem  in  Chapter  V  and  some  of  the  numerical  results  are 


listed. 
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CHAPTER  II  MEASURE  OF  ACCURACY 

2.1  INTRODUCTION 

As  stated  in  section  1.4,  we  are  investigating  limitations  on  the 
accuracy  of  stiff  methods  (definition  in  section  1.2).  The  results  of  the 
investigation  are  meaningful  only  if  a  reasonable  measure  of  accuracy  is 
used.  For  all  practical  purposes,  we  would  say  that  numerical  method  A  is 
more  accurate  than  method  B  if  the  numerical  results  from  using  method  A 
are  closer  approximations  to  the  true  solution  than  those  from  using  method 
B  for  the  same  stepsize.  Consequently,  the  order  of  the  method  alone  is 
too  crude  as  a  measure  of  accuracy;  rather  a  more  refined  measure  of  global 
error  is  needed.  Since  the  actual  global  error  accumulated  will  depend  on 
the  differential  equation  to  be  solved,  parameters  which  depend  only  on  the 
numerical  method  in  either  a  global  error  bound  or  in  an  expression  for  the 
asymptotic  global  error  will  be  suitable  as  a  measure  of  how  accurate  the 
method  is. 

Two  candidates  for  the  measure  of  accuracy  are  introduced  in  sec¬ 
tion  2.2  and  2.3,  respectively.  To  simplify  expressions  for  global  error, 
we  assume,  for  the  rest  of  this  chapter,  that  the  starting  values  are  ex¬ 
act,  and  that  the  roundoff  errors  are  negligible  compared  with  the  discre¬ 
tization  errors. 

2.2  ASYMPTOTIC  GLOBAL  ERROR 

We  first  consider  a  measure  of  accuracy  which  indicates  the  mag¬ 
nitude  of  the  asymptotic  global  error,  i.e.,  the  magnitude  of  the  most  dom¬ 
inant  term  as  the  stepsize  h  becomes  arbitrarily  small  (Gear,  1971, 
pp. 75-6, 204-5): 
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yn  -  y(tn)  =  hkck+1  ftn  l<(tn,t)y(k+1)(T)dT  +  0(hk+1)  h  -  0  (2.1) 

J  t0 

where  K(t,x)  satisfies  the  homogeneous  ordinary  differential  equation 
3  K  3f 

— (t.x)  =  — (t.y(t))K(tir)  K(t,t)  =  I 

at  3y 

Thus,  one  measure  of  accuracy  is  the  parameter  ck+-|,  which  is  the  only 
method  dependent  part  of  the  leading  term  of  the  asymptotic  error  expan¬ 
sion.  The  error  constant  can  be  expressed  in  terms  of  the  a.'s  and  the 

J 

B.'s  as 
J 

Ck+1  -  — —  S  [jk+1oi  -  (k-H)j'V]  (2.2) 

(k+l)!  j=l  3  3 

and  in  terms  of  the  b.'s  as  (Genin,  1973) 

J 


1 

ck+l  =  " 

j 


j=0  j+1 
even 


(2.3) 


2.3  GLOBAL  ERROR  BOUND 

Under  the  assumption  that  the  stepsize  h  used  when  applying  (p, a) 
to  solve  system  (1.1)  is  small  enough  such  that 


Henrici  (1962,  p.248)  shows  that  an  upper  bound  for  the  global  error  is 
given  by 
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(tn-to)rBL 

k  (tn't0)eXPC  l-h|y0|L  ] 

iy»-y(t»)lihrGY - i-m^|L-P . 


(2.4) 


where 

?0  =  6k/c,k 

|y(k+1)(t)|  <  Y  ts[t0,T] 

1  “  i 

r  =  max |y •  |  ,  - r - rri - •=  s  (2.5) 

j>_0  J  aQ<;  +a-|C  +...+ak  j=0  J 

k 

B=  I  1 3  - 1 
j=0  J 

rk 

G  =  (  1 G ( s ) I ds  (2.6) 

J  0 

G(s)  is  the  influence  function  associated  with  (p,a)  and  is  defined  by 

G(s)  =  —  E  [Wj-s)k  -  kfUj-s)*'1]  (2.7) 

k!  j=0  3  3 

where,  for  any  ve(-ooj00)t 


x  >  0 
x  <  0 


Mote  that  a  method  of  order  >_  k  is  exact  for  polynomials  of  degree  k  and 
hence  G(s)  =  0  for  s  <  0.  Clearly,  by  definition,  G(s)  =  0  for  s  >  k.  It 
is  easily  shown  using  (1.6)  that  G(s)  is  the  kernel  of  the  operator  l_h 
(Henrici,  1962,  p.248): 


Lhy(t)  =  hk+1  f  G(s)y^k+1^(t+sh)ds 

•'—co 


(2.8) 


where  Lh  is  the  linear  difference  operator  defined  by 
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k  k 

Lhy(t)  =  e  ot^y (t+ih)  -  h  e  BiY*  (t+ih)  (2.9) 

n  i=0  i=0 

Henrici  (1962,  p. 242)  proves  that  r  is  finite. 

A  refinement  of  the  above  bound  can  be  derived  by  using  the  quan¬ 
tity 


e  rk+B  rk-T+.-.+e 


r  =  max  |  y-s  I  • 
j>0  J 


=  E  Y.C- 


a0?k+al?k~^+“,+ak  J=0  ^ 


(2.10) 


instead  of  rB,  as  follows: 

We  start  with  the  difference  equation  satisfied  by  the  numerical 

solution, 

k  k 

Vj-kM  -  h  ■  0  J  i  k 

The  difference  of  the  above  equation  and  (2.9),  the  latter  evaluated  at  t  = 
t._k,  yields 


•iej-k+i  -  h  eisj-k+i  =  -  Lh^(tj-k> 


j  >  k 


where  e^  =  yy  -  y(ty)  and  §v  =  f(tv,yv)  -  f(tv,y(tv))  for  all  v  >  0.  In 

particular,  e^  =  e^  =  0  for  0  1  v  £  k-1.  Multiplying  the  above  equation  by 

y  .  and  summing  from  j  =  k  to  j  =  n, 
n  j 


“iej-k+i 


A  Yn'j  i=o  ei’eJ'k+1 


The  first  term  is,  by  (2.5),  simply  e  ,  whereas  the  second  term  can  be 
shown,  using  (2.10),  to  be 
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h  A  Wj 


n-k 


Hence,  we  have 

en  '  h  jk  Vj6j  '  jf0  Y0Lhy(t"-k-j> 

As  in  Henrici  (1962,  p.247),  Lhy(t)  can  be  uniformly  bounded  by 

|Lhy(t:)|  <  hk+1GY  tQ<t<t+kh<T 

From  (2.11),  by  the  triangle  inequality, 

n-1 

lenl  1  h?L  .1  lejl  +  h|y0IMen|  +  hk  rGY(n-k+l ) 
j=k 


(2.11) 


or 


h?L  n-1 

Ienl  <  -  Z 

n  - l-h|y0|L  j-k 


h  rGY 

|  e .  |  + - (n-k+1 ) 

J  l-h|^olL 


if  h | yq | L  <  1. 

We  shall  make  use  of  the  discrete  Bellman  inequality  as  stated  in 
the  following  lemma  (Babuska,  Prager,  Vitasek,  1966,  p.57): 


Lemma  2.1 


Let  <t  ,  ib  ,e  be  sequences  defined  for  v  =  0,1,..., n,  0  >  0, 

Yv  v  v  v  — 

and 

v-1 

^  ^  +  z  e  ^  v  =  0,1,. ..,n 

v  -  v  y_Q  y  y 


Then 
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n-1  n-1 

+  *  Vm  ”  <1+®s> 

y=U  S=y+1 

Applying  the  above  lemma  with  <j>v  =  |  I »  and 


k+1 

h  rGY  hn 

ip  = - v,  9  = - 

v  1  “h | y0  I  L  V  l-h|y0|L 


we  obtain  the  following  global  error  bound: 

expc  — — — ]  -  1 

k  1-MyqIl 

|yn  -  y(tn)l  1  h  rGY - r -  (2.12) 

rL 

The  bound  (2.12)  is  better  than  Henrici's  bound  (2.4)  since 

f  <  rB 

The  inequality  is  usually  strict.  For  example,  for  the  2-step  Adams  Bash- 
forth  method  (Gear,  1971,  p. 1 09) , 

f  =  3/2 
rB  =  2 

From  the  bound  (2.12),  one  possible  measure  of  accuracy  is  the 
quantity  rG.  There  is  no  nice  expression  for  rG  in  terms  of  the  b.'s. 

J 

We  remark  that  the  bound  (2.12)  can  be  further  improved  by  using 

the  fact  that,  from  (2.8), 

n-k  k+1  ±  n-k 

|  z  Y-;  Lhy(t  k_.)l  lh  Y  (  I  s  YiG(s+j)|ds 
j=0  J  n  n  K  J  Jk-n  j=0  J 


If  we  define 


1  rk  n-k 

=  max  - 1  |  z  Y1-G(s+j)|ds 

k£n£N  n-k+1  J  k-n  j=0  J 


x 
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then  using  a  similar  argument  as  in  deriving  (2.12)  from  (2.11),  we  obtain 
a  global  error  bound  similar  to  (2.12)  except  that  the  quantity  rG  is  re¬ 
placed  by  x.  The  fact  that  x  1  rG  is  obvious,  and  x  should  be  an  impro¬ 
vement  for  methods  having  an  influence  function  G(s)  which  changes  sign  and 
y.'s  approximately  the  same. 

J 

2.4  CONCLUDING  REMARKS 

We  have  introduced  two  possible  candidates  as  reasonable  measures 
of  accuracy,  namely,  | c k+-|  |  and  rG  (or  x).  It  is  not  difficult  to  show 
that 


Moreover,  equality  holds  if  G(s)  is  of  the  same  sign  on  (0 , k) . 

The  asymptotic  error  (2.1)  could  be  too  optimistic  because  of  its 
asymptotic  nature,  whereas  the  error  bound  (2.12)  could  as  well  be  too  pes¬ 


simistic.  In  fact,  ck+-|  depends  only  on  the  even  bj's,  indicating  possible 
inadequacy  as  a  useful  measure  of  accuracy.  On  the  other  hand,  rG  (or  x) 


is  not  the  only  method  dependent  quantity  in  the  error  bound,  and  there  is 
one  more  important  quantity  in  (2.12),  namely,  r.  Because  |ck+-|  |  can  be  ni¬ 
cely  expressed  as  a  linear  combination  of  the  b.'s  (cf.  (2.3)),  we  shall 

J 

choose  it  as  the  measure  of  accuracy  in  our  study. 


It  should  be  finally  remarked  that  the  measure  |ck+-||  has  the  de¬ 
fect  that  it  does  not  reflect  the  instability  of  the  method  (p,p),  in  the 


sense  that  its  magnitude  could  be  small  even  though  the  method  is  not  sta¬ 
ble  at  all.  Hence  it  is  a  measure  of  accuracy  for  stable  methods  only.  It 


would  be  desirable  to  have  a  measure  of  accuracy  that  also  includes  the  ef¬ 
fects  of  the  positions  of  the  roots  of  p(?),  since  where  the  roots  lie 
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could  have  an  effect  on  the  numerical  solution,  as  in  the  case  when  there 
are  multiple  roots  near  the  unit  circle.  If  we  consider  the  case  when  the 
roots  of  p(?)  are  real  and  have  magnitude  greater  than  or  equal  to  1-e, 
ee(0,l),  then  r  will  be  large  if  e  is  small.  In  fact, 

j  k+v-2  1 

r  >  max  z  (  )(l-e)v  =  "77 

j>0  v=0  k-2  e  ~ 

Hence,  the  measure  rG  (or  x)  is  preferable  in  that  sense. 
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CHAPTER  III  FROM  Ag-STABILITY  TO  A(a)-STABILITY 
3.1  INTRODUCTION 

As  pointed  out  in  section  1.4,  the  reason  for  considering  Ag- 
stability  and  A(a)-stability,  which  are  weaker  than  A-stability,  is  to  al¬ 
low  for  methods  of  order  greater  than  2.  In  proving  the  existence  of 
Ag-stable  methods  in  L[ k]  (definition  in  section  1.3)  for  arbitrary  k, 
Cryer  (1973)  shows  that  the  following  methods,  as  parameterized  by  the  po¬ 
lynomial  s(z)  (cf.  section  1.3), 

s(z)  =  (z  +  d)k 


are  Ag-stable  for  d  >_  2k+1.  Jeltsch  (1976)  later  proves  that  these  methods 
are  in  fact  A(0)-stable.  Hence  A(0)-stable  k-step  methods  of  order  k  exist 
for  arbitrary  k.  The  error  constant  of  Cryer's  method  is  given  by  (cf. 
(2.3)) 


1  k  k  dk-j 

—  z  (  ) - < 

2K  j=0  j  j+1 
j  even 


I  C- 

and  thus  |ck+1|  is  greater  than  2K  .  Although,  as  remarked  by  Cryer,  the 
bound  for  d  could  be  strengthened,  yet  the  essential  point  is  that  ck+-|  is 
of  order  0(dk),  and  that  as  k  increases  so  must  d.  It  is  the  impractica¬ 
bility  of  Cryer's  methods  that  motivates  our  investigation  of  limitations 
on  the  accuracy  of  Ag-stable  methods. 

Some  characterizations  and  properties  of  Ag-stable  methods  are 
listed  in  section  3.2.  In  section  3.3,  we  discuss  in  detail  the  accuracy 
of  Ag-stable  methods  in  l[ 1]  and  L[ 2].  Some  lower  bounds  on  |ck+1|  for 


Ag-stable  methods  in  L[k],  for  arbitrary  k,  are  derived  in  section  3.4. 

The  existence  of  A(a)-stable  methods  of  arbitrary  order  k  for  any  cxe(0,Tr/2) 
is  established  in  section  3.5. 

3.2  CHARACTERIZATION  AND  PROPERTIES  OF  AQ-STABLE  METHODS 

This  section  contains  a  survey  of  known  results  concerning  AQ- 
stability  which  will  be  referred  to  in  subsequent  sections  and  also  in 
Chapter  V.  The  parameterization  (r,s)  for  methods  in  L[k]  (cf.  section 
1.3)  is  used  throughout. 

We  characterize  AQ-stable  methods  in  terms  of  some  algebraic  con¬ 
ditions  on  the  polynomials  r(z)  and  s(z).  Restrictions  upon  the  coeffi¬ 
cients  b.'s  are  also  given.  All  methods  are  assumed  to  be  in  L[ k].  For 
any  complex  number  z,  Arg  z  denotes  that  value  of  arg  z  which  is  in  the  in¬ 
terval 

We  first  consider  AQ-stabil ity.  The  following  theorems  are  taken 
from  Cryer  (1973): 

Theorem  3.1 

The  following  statements  are  equivalent: 

(i)  (r,s)  is  AQ-stable. 

(ii)  For  all  qe(-~,0),  the  roots  of  r(z)-qs(z)  lie  in  the  interior  of  the 
left  half  complex  plane. 

(iii)  For  Re  z  >  0,  r(z)/s(z)  is  regular,  and  does  not  take  values  in 
(-oo,0).  For  z  =  iw,  i  -  AT,  we(-“>,“),  r(iw)s(-iw)  does  not  take 
values  in  (-°°,0). 

Theorem  3.2 

If  (r,s)  is  Ag-stable,  then  the  zeros  of  r(z)  and  s(z)  lie  in  the 
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closed  left  half  plane,  and  those  on  the  imaginary  axis  are  at  most  double 
zeros. 

Theorem  3.3 

If  (r,s)  is  Ag-stable,  then 

b0  1  0,  b-]  1  0 


and,  if  k  i  3, 


bj  >  0  for  j  =  2,3,...,k-l 

A  characterization  of  A(a)-stable  methods  is  given  by  Widlund 


(1967): 


Theorem  3.4 

The  following  statements  are  equivalent: 

(i)  (r,s)  is  A(a)-stable. 

(ii)  For  all  qe$a  =  (z  :  |Arg(-z)|  <  a,  z  4  0>,  the  roots  of  r(z)-qs(z)  lie 
in  the  interior  of  the  left  half  plane. 

(iii)  For  Re  z  >  0,  r(z)/s(z)  is  regular,  and  takes  its  values  in  the  com¬ 
plement  of  S  . 

a 

If,  in  addition  to  A(a)-stability,  we  assume  that  (r,s)  is  A^- 
stable,  then  r(z)/s(z)  is  regular  on  Re  z  0.  Thus,  by  continuity, 
A(a)-stabil ity  implies  that 

r(iw)/s(iw)  ^  for  all  we(-°°,°°)  (3.1) 

On  the  other  hand,  the  locus  r(iw)/s(iw),  we(-°°,°°),  of  an  A^-stable  method 
divides  the  complex  q-plane  into  several  components,  one  of  which  contains 
a  neighborhood  of  °°,  denoted  by  V.  Obviously,  if  (3.1)  is  satisfied,  V 
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contains  Sa.  Since  r(iw)/s(iw),  we{-®,“),  is  the  locus  of  q  in  the  complex 
plane  for  which  a  root  of  r(z)-qs(z)  is  purely  imaginary,  all  roots  of 
r(z)-qs(z)  for  any  q  in  the  interior  of  V  have  to  lie  in  the  interior  of 
the  left  half  complex  plane,  a  condition  satisfied  by  r(z)-qs(z)  at  q  =  ». 
From  theorem  3.4,  this  means  that  the  method  is  A(a)-stable.  Hence,  for 
methods  which  are  Aro-stable,  (3.1)  is  necessary  and  sufficient  for  A(a)- 
stability. 

It  can  easily  be  shown  that  (3.1)  is  equivalent  to  the  following 

condition: 


Re  r(iw)/s(iw) 

-  - £  cot  a  (3.2) 

| Im  r(iw)/s(iw)| 

for  all  we(-oo,o°)  such  that  r(iw)  4  0*  Condition  (3.1)  or  (3.2)  will  be 
used  in  subsequent  sections  and  in  Chapter  V  to  prove  A(<*)-stability  when 
the  method  is  A^-stable. 

We  conclude  this  section  by  giving  another  necessary  and  suffi¬ 
cient  condition  for  A(a)-stability  if  the  method  is  known  to  be  Ap-stable 
instead  of  Aro-stable.  The  condition  is  a  direct  consequence  of  a  theorem 
by  Jeltsch  (1976),  rewritten  using  the  parameterization  (r,s)  instead  of 
(p,a): 


Theorem  3.5 

The  conditions  (i)-(iv)  are  necessary  and  sufficient  for  a  method 
in  L[k]  to  be  A(0)-stable: 

(i)  (r,s)  is  Ag-stable. 

(ii)  The  roots  of  s(z)  on  the  imaginary  axis  are  simple. 

(iii)  Let  z  =  iw  be  a  purely  imaginary  root  of  r(z),  then 
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s(iw) 

Re - >  0 

r'(iw) 

( i v )  Let  z  =  iw  be  a  purely  imaginary  root  of  s(z),  then 

r(iw) 

Re  - >  0 

s'(iw) 

If  the  method  (r,s)  is  Ag-stable,  then  it  is  A(a)-stable  if  and 
only  if  conditions  (ii)-(iv)  in  the  above  theorem  hold,  and  (3.2)  holds  for 
all  we(-c»,oo)  such  that  r(iw)  f  0. 

3.3  SPECIAL  CASES  k  =  1,2 

The  limitations  of  the  accuracy  of  AQ-stable  methods  in  £[1]  and 
£[ 2]  are  discussed  in  this  section.  A  thorough  analysis  is  possible  since 
only  a  few  parameters  are  involved. 

3.3.1  k  =  1 

The  polynomials  r(z),  s(z),  and  the  error  constant  c2  of  a  method 
inL[l],  parameterized  by  (bQ),  are 

r(z)  =  2 
s(z)  =  z  +  b0 

c2  *  -  V2 

Since 

r(iw)  2bn  2w 

- =  —z — 2  "  — o  for  all  we(-“,“») 

s(iw)  w  +bj  w  +bg 

the  method  (bQ)  is  A-stable  if  and  only  if  bQ  0.  If  bQ  =  0,  we  have  the 
trapezoidal  rule,  which  is  in  fact  of  order  2  since  c2  =  0.  Its  region  of 
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And  hence 


s  (z )  =  z2  +  b-jZ  +  b( 


1 


1 


-  -  -(Oq  +  ~ ) 


4 


r(iw)  b^Q  -  iw(w2+b2-bQ) 
s ( i w )  (bg-w2)2  +  b2w2 


(3.3) 


By  theorem  3.3,  if  (r,s)  is  Ag-stable,  it  is  necessary  that 


b0  >  0,  b1  >  0 


Thus  there  exists  no  Ag-stable  method  in  l[ 2]  with  |c3|  <  1/12.  The  po¬ 
lynomial 


r(z)  +  qs(z)  =  qz2  +  (qb1+2)z  +  (qbQ+2b1 ) 
has  roots  in  the  open  left  half  plane  for  any  qe(0,°°)  if  bg,^  satisfy 

bg  >_  0,  b-j  >_  0,  bg+b-j  >0  (3.4) 

Thus,  from  theorem  3.1,  any  method  (bg,^)  which  satisfies  (3.4)  is  Ag- 
stable.  If  bg  =  0,  b-j  >_  0,  r(z)  and  s(z)  have  a  common  factor  and  the 
corresponding  method  reduces  to  the  one-step  trapezoidal  rule.  Note  that 
it  is  the  only  Ag-stable  method  whose  order  (2)  is  greater  than  its  step 
number  (1)  (Cryer,  1973). 


r(iw) 

Re - 

s(iw) 


2bl  bg 


/l  -  2 \ 2 . .  2  2 

(bg-w  )  +b-j w 


From  (3.3), 
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is  nonnegative  for  any  we{-~,“)  if  bg  >_  0,  and  b-j  >_  0.  Hence  for  any  |Cg| 
>  1/12,  we  can  construct  A-stable  methods  in  z,[2]  by  choosing  bg  -  4 1 c^l  - 
1/3,  b i  ^  0  (cf.  (3.2)  and  the  remarks  after  theorem  3.5).  If  by  >  0,  the 
method  is  A  -stable. 

oo 

The  polynomials  p(c)  and  o(?)  corresponding  to  (bg.b-j)  are 


b,+l  0  b,-1  b,+l 

— — s2  -  b,c  ♦  -1—  ■  -L-(c 
2  1  2  2 

yv^2  +  1-bo  +  bo~bi+1 


b-i-l 

DU  -  — ) 


so  that 


y .  =  1  -  (— )J+1  for  all  j  >  0 

J  bi+1 


l-bl 

-s(s 

+ 

1 

0  <  s  <  1 

4 

V1 

1+bl 

(s  - 

2)(S  -  1  + 

-Q~) 

1  <  s  <  2 

4 

1+b, 

0 

otherwise 

For  fixed  >_  0,  if  we  choose  0  £  b^  £  1+bg,  then  ®(s)  1S  one  S19n  on 
[0,2]  and  hence  G  =  |c3| ,  and  if  we  choose  1  <  b^  then  r  £  1.  In  either 
case,  the  method  is  A-stable,  and  is  also  A^-stable  if  bQ  >  0  and  b]  >  0. 

3.4  SOME  LOWER  BOUNDS 

As  pointed  out  in  section  2.2,- the  error  constant  ck+1  associated 
with  a  method  in  L[k],  parameterized  by  ( bg , b-j  ,...,bk_-j ),  can  be  expressed 
linear  combination  of  the  even  b  j ' s  (cf.  (2.3)): 


as  a 
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'k+1 


1  k  b. 


2  j=0  j+1 

j  even 


In  this  section,  we  derive  some  lower  bound  on  jc^il  for  methods  inl[k] 
which  are  Ag-stable. 

We  first  consider  the  case  when  k  is  even.  A  necessary  condition 
for  the  method  (bg,b-j,...,b|<_i)  to  be  Ag-stable  is  that 


bg  1.  0» 


b1  >  0 


bj  >  0  j  =  2,3,...,k-l 


(cf.  theorem  3.3).  Since,  by  normalization,  bk  =  1,  a  lower  bound  for 
Ick+il  is 


'W  ■■'wSVii 

the  inequality  being  strict  if  k  >_  4.  The  case  when  k  =  2  has  been  dis¬ 
cussed  in  section  3.3.2,  where  it  was  shown  that  the  lower  bound  is  at¬ 
tained  by  the  trapezoidal  rule. 

A  similar  argument  when  applied  to  the  case  when  k  is  odd  will 
lead  to  the  lower  bound 

I ck+l I  =  "  ck+l  -  0 

with  strict  inequality  when  k  3.  Again,  the  lower  bound  is  attained  by 
the  trapezoidal  rule  in  the  case  when  k  =  1  (cf.  section  3.3.1).  It  will 
be  shown  in  the  next  section  that  for  any  c4  <  0,  there  exist  methods  in 
Z/[3]  which  are  Ag-stable  and  have  error  constant  c4.  For  the  general  case 
when  k  1  5,  k  odd,  a  lower  bound  for  |ck+-|l  can  be  derived  from  the  fol¬ 
lowing  lemma: 
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Lemma  3.1 

Let  k  >5,  k  odd.  If  (b()tb1#...tbk_1)  is  an  A0-stable  method  in 

L[ k],  then 


r 


where 

Hk(q)  =  1 

H  •  (q )  =  b.  +  q—  j  =  0,1,...,  k-1 
J  J  2 

lie  in  the  interior  of  the  left  half  complex  plane.  Hence,  it  is  necessary 
that  the  following  inequality  be  satisfied  for  all  q  >0  (cf.  Appendix  A): 

Hk_i(q)H-j(q)  -  Hg(q)  >  0 


In  terms  of  powers  of  q,  the  left  hand  side  can  be  written  as 
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a.  0  a..^  ..  an 

7" +  (r_  bi)]q  +  (bk-ibi  •  bo> 

Since  b  >  0  for  v  :>  2  (cf.  theorem  3.3),  the  coefficient  of  q  is  negative 
from  (3.5).  In  order  that  the  above  quadratic  (in  q)  be  positive  for  all  q 
>  0,  it  is  necessary  that 

-  (7-  bi)]2  *  2ai  (bk-ibi  •  V 


2a,b0  +  -  -)2  -  2(i^l  -  A  -  b,  -  -)  + 

u  2  k  2  k  2  k 


+  ^  "  bl  "  7?2  ‘  2albk-lbl  <  0 


(3.6) 


which  implies  that 

<7'bi  - ->2  -  2a1bk-lb!  <0 

since  the  first  three  terms  in  (3.6)  are  nonnegative.  However,  from  (1.9), 


\2  bk-2b3 


/  U  l.  \C  Is- 1.  O  • 

( - b,  -  >  -  >  - D,  > 

2  1  k  “  3 ( k-2 )  “  1 2 ( k-2 )  1  “ 


k  =  5 


k  >  7 


i  2aibk-ibi 

making  use  of  the  fact  that,  by  theorem  3.2,  s(z)  has  no  root  in  the  open 
right  half  plane,  and  hence  (cf.  Appendix  A) 

4bk.2b3  >  (k-1 )2b, 


We  thus  have  a  contradiction. 


Q.  E.  D. 
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if  the  method  is  AQ-stable.  Moreover,  the  infimum  of  |c^+^|  for  AQ-stable 

2 

methods  lies  between  l/(2k3k)  and  2k  +k~^. 

3.5  HIGHER  ORDER  A(a)-STABLE  METHODS 

In  this  section,  we  demonstrate  how,  for  given  ae  (0,ir/2) , 
A(a)-stable  methods  in  l[ k]  can  be  constructed,  and  by  doing  so,  prove  that 
A(a)-stable  methods  of  arbitrary  order  k  exist  for  any  ae(0,Tr/2).  The  ac¬ 
curacy  of  such  methods  is  also  discussed.  A  statement  about  the  accuracy 
of  AQ-stable  methods  ini [3]  is  given  as  a  corollary.  Since  the  cases  when 
k  =  1 ,2  have  been  analyzed  in  section  3.3,  we  assume  throughout  that  k  >_  3. 

Dill  and  Gear  (1971),  when  searching  for  stiffly  stable  methods 
of  order  greater  than  6  (cf.  section  1.4),  make  an  interesting  observation 
that  there  is  a  greater  likelihood  of  stiff  stability  for  methods  whose 
p(c)  polynomial  has  multiple  roots  near  1.  Consider  the  method  represented 
by 

P(0  =  (c-l)k,  a(c)  «  (c+6)(c-Dk_1/(l+6) 

The  common  factor  (c-l)k_1  makes  the  method  of  order  k  automatically.  The 
method  "almost"  satisfies  the  root  condition  and  has  the  same  stability  re¬ 
gion  as 

pfe)  -  c-l.  0(c)  ■  (c  +  $)/(l  +  6) 

which  is  A-stable  if  6e(-l,l].  Now  by  choosing  k-1  roots  of  a(c)  to  be 
slightly  less  than  1  instead  of  1 ,  we  may  be  able  to  get  an  A(a)-stable 
method  in  l{ k]  with  a  close  to  tt/2.  To  this  end,  consider  the  two 
parameter  family  of  methods  parameterized  by  the  following  polynomial: 

k-1 


s  (z )  =  (z  +  d)(z  +  D)1 


(3.7) 
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U-, (z)  =  (u(z)  -  u(0))/z 

Since  s-j(z)  =  u(z)  +  du-j(z),  the  rational  function  r(z)/s(z)  for  z  =  iw, 
is  then  given  by 

r(iw)  2  Ui(iw)  k  sm(iw) 

— —  - D+d - +  s  - ]  (3.12) 

s(iw)  iw+d  u(iw)  m=3  mu(iw) 

m  odd 

If  we  are  able  to  choose  d  and  D  so  that  the  quantity 

u i ( i w )  k  s  (iw) 

d - +  x  - 

u(iw)  m=3  mu(iw) 
m  odd 

is  sufficiently  small,  for  any  we(-°°,<»),  with  respect  to  the  constant  1, 
then 


r'(iw)  _  2 

- =  -  for  all  we (-”,“) 

s(iw)  iw+d 

in  which  case  we  may  be  able  to  get  A(cc)-stability  for  a  close  to  tt/2.  And 
that  is  the  main  idea  behind  the  discussion  that  follows.  We  shall  find  a 
uniform  upper  bound  in  terms  of  d  and  D  for  the  following  "undesired"  por¬ 
tion  of  the  function  r(iw)/s(iw)  over  -»  <  w  =  Dy  <  °°: 

Ui(iDy)  k  s„(iDy) 

d - +  ^  - 

u(iDy)  m=3  mu(iDy) 
m  odd 

Then  by  appropriately  choosing  d  and  D,  it  will  be  shown  that  we  are  able 
to  get  A(a)-stabil ity  for  any  given  ae(0,ir/2).  Note  that  2/(iw+d), 
we(-co,oo)j  is  the  locus  of  a  circle  in  the  complex  plane  centered  at  1/d 
with  radius  1/d. 

An  expression  for  | u ( i Dy )  |  is  easily  obtained  from  (3.7)  and 
(3.11)  by  substituting  z  =  iDy  and  factoring  out  the  factor  D, 
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k-1 

| u ( i Dy ) |  =  Dk"^(y2  +  1)  2  (3.13) 

We  next  find  an  upper  bound  for  | sm(i Dy ) | ,  0  £  m  <  k.  From 
(3.8),  (3.9),  and  (3.10), 

Un,(TDy)lz  *  Ik"{  Dk"n  'iBnH.j(fDy)Jl2  i  D2(k"n)pk.n,(y2) 

J=0 

2  2 
where  Pk_m(y  )  is  a  monic  polynomial  of  degree  k-m  in  y  ,  independent  of 

the  parameters  d  and  D.  Combining  (3.13)  and  the  above  bound,  we  have,  for 

any  m  >_  1 , 

,sm(iD^|  < 

u(1Dy)  “  D"1'1 

uniformly  in  ye(-°°,°°),  in  which  M'  is  some  positive  constant.  It  follows 
that,  if  D  1 , 

M 


k  s  ( i Dy )  M'  k  1 

^  111 _  <  _  2  _ _  <  — 

m=3  mu(iDy)  D2  m=3  mDm  ^  D2 


(3.14) 


m  odd  m  odd 

where  M  is  a  constant  independent  of  d  and  D  (e.g.  M  =  M'k). 

By  a  similar  argument,  we  can  prove  that  there  exists  a  constant 
K  >  0  such  that 

u-i(iDy)  K 

l~J - 1  1  - 

u(iDy)  D 

uniformly  in  ye(-»,»).  Together  with  (3.14),  we  obtain  the  following  bound 
far  the  "undesired"  portion  of  r(iDy)/s(iDy)  (cf.  (3.12))  when  D  >_  1: 


Ui(IDy)  k  s  (iDy)  dDK  +  M 

|  d - +  Z  - 1  1 - 2 — 

u(iDy)  m=3  mu(iDy)  D 

m  odd 


(3.15) 


We  now  proceed  to  show  how,  for  any  given  ote(0,7r/2),  the 


40 


parameters  d  and  D  can  be  chosen  so  that  the  method  is  A(a)-stable.  Assume 
that  ae(0,7r/2)  be  given.  To  get  A(a)-stability,  we  want 

r(iDy) 

I  Arg  - 1  <  ir  -  a 

s(1Dy)  " 

for  any  ye(-°°,oo)  such  that  r(iDy)  =f  0  (cf.  (3.1)).  Consider  an  arbitrary 
ye(-«>foo).  if  we  choose  d  and  D  >_  1  such  that 

dDK  +  M 

- 5 — £  sin(-  -  a)  (3.16) 

2 

then 


r(iDy)  2 

I  Arg  - 1  <  |  Arg - 1 

s(iDy)  iDy+d 


Ui(iDy)  k  s  (iDy) 

+  | Arg  [1  +  d-! - +  z  - ]| 

u(iDy)  m=3  mu(iDy) 
m  odd 


-  +  sin-1  (- 
2 


dDK  +  M 


-)  <_  ir  -  a 


Therefore,  from  (3.1),  the  method  (r,s),  if  convergent,  is  A(a)-stable  pro¬ 
vided  (3.16)  holds.  It  remains  to  show  that  the  polynomial  r(z)  ,  when 
(3.16)  is  satisfied,  has  roots  in  the  left  half  plane.  We  shall  need  the 
following  theorem  (Levinson  and  Redheffer,  1970,  p. 21 8 ) : 


Theorem  3.7  (Rouch^'s  theorem) 

Let  f(0  and  g(c)  be  analytic  in  a  simply  connected  domain  con¬ 
taining  a  Jordan  contour  J.  Let  |f(c)|  >  |g(c)l  on  J.  Then  f(?)  and 
f(s)+g(?)  have  the  same  number  of  zeros  inside  J. 


r(z)  =  2[u(z)  +  du-j 


(2)  + 


k 

z 

m=3 


m  odd 


Since 
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from  (3.15)  and  (3.16), 

r(iw)  , 

| - u(iw)|  <  |u(iw)|  (3.17) 

for  any  we(-®,°°).  Applying  theorem  3.7  with 

C+l  1  s+1  s+1 

f(c)  =  U(— ).  gU)  =  -r( — )  -  u (  ) 

S-l  2  s-1  c-1 

and  J  being  the  unit  circle,  we  conclude,  from  (3.17),  that  the  polynomials 

u(z),  r(z)/2 

have  the  same  number  of  zeros  in  the  interior  of  the  left  half  z-plane. 
Hence  all  the  roots  of  r(z)  lie  in  the  open  left  half  plane,  and  the 
corresponding  method  (r,s)  is  convergent  (cf.  section  1.3)  —  in  fact, 
strongly  stable. 

We  can  hence  obtain,  for  any  ae(0,fr/2),  methods  in  L[k]  which  are 
strongly  stable,  A^-stable,  and  A(a)-stable  by  choosing  the  parameters  d 
and  D  >.  1  such  that  (3.16)  is  satisfied  (Slight  modification  needed  if  D  < 
1).  It  should  be  remarked  that  d  can  be  made  arbitrarily  small,  and  in 
fact  zero,  though  the  method  will  not  be  A^-stable  when  d  =  0. 

From  (3.8),  the  magnitude  of  the  error  constant  ck+-j  is  given  by 


k  k-1  Dk'J 

e  (  ) - +  d 

j=2  j-1  j+1 
even  j 


k-1 

z 

j=0 


k-1  Dk'J_1 

(  ) -  ] 

j  0+1 


even 


(3.18) 


If  we  choose  d  such  that  d  =  0(D"  ),  then  from  (3.16),  for  A(a)-stability, 

D2  =  0(- - )  as  a  +  2 


Hence,  the  magnitude  of  the  error  constant,  being  of  order  0(D  )  for  lar- 
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ge  D,  is  of  order 

,  H 

0([- - ]  2  )  as  «  -  I 

|  -  a  2 

The  results  in  this  section  can  be  summarized  by  the  following 

theorem: 

Theorem  3.8 

For  arbitrary  k  and  any  ae(0,7r/2),  there  exist  A(a)-stable  k-step 
methods  of  order  k  which  are  AM-stable  and  strongly  stable. 

Corollary  3.8.1 

The  method  corresponding  to 

s(z)  =  z(z  +  D)k~^  (3.19) 

is  Ag-stable  if  D2  l  2k  3  +  1. 

Proof 

To  prove  the  corollary,  all  we  have  to  do  is  to  obtain  value  for 
M'  (cf.  (3.14))  in  terms  of  k  when  k  >_  3  (the  case  when  k  <  2  is  trivial  -- 
see  section  3.3).  We  shall  make  use  of  the  following  identities,  the  veri¬ 
fication  of  which  is  straightforward: 

k-1  k-1 

^  ( 

j=0  j 

j=0  mod  4 

k-1  k-1 

^  ( 

j=l  J 

j=l  mod  4 


k-3  /, 

k-3  2  'u 

)  =  2k  13  +  2  c  cos - 


(3.20a) 


k-3  n 
k-3  ~T 

)  =  2K  J  +  2  ^  sin - 


4 


(3.20b) 
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j=3  j 
j=3  mod  4 

Since  for  any  0  <_j  <_  k-1 , 


k-1  k-1  i<_3  2  (k-1)  it 

E  (  )  =  2K  -  2  cos - 

i=2  j  4 

j=2  mod  4 

k-1  k-1  n  y  (k-1)  ir 

z  (  )  =  2  -2  sin - - 


(3.20c) 


(3.20d) 


|y|J  i  (y  +  i) 


for  all  ye(-oo,co) 


we  have,  for  m  >  3,  m  odd,  from  (3.8),  (3.10),  (3.13),  and  (3.20), 

,  x  s  (iDvl  „  k-m  k-1  k-m  k-1 

,  r»,.  r  (  ).  E  (  )}] 


(iDy)  o  k-m  k-1  k-m  k-1  « 

m  |2<rmaxf  z  (  )•  S  (  )>]  + 

u(iDy)  ”  3=0  m+j-1  3=2  m+j-1 

j=0  mod  4  3=2  mod  4 


k-m  k-1  k-m  k-1  . 

+  [max{  E  (  )»  ^ 

j*l  m+j-1  3=3  m+j-1 

j=l  mod  4  j=3  mod  4 


<  2[2k-3  +  2  2  ]2  <  [2  2]2 


Therefore,  if  D  _>  1, 


k-1 

2  2  2k"3 


k  s  (iDy)  2  k  1  2  2 

|  £  - 1  <  — 5-  z  — £ - o -  <  “5 — 

m=3mu(iDy)  D2  m=3  mD"1-3  3(D2-1)  D  -1 

m  odd  m  odd 

uniformly  in  ye{-»,»).  We  note  that  the  above  upper  bound  is  different 

k_  3/2 

from  that  in  (3.14).  In  fact  (3.14)  is  satisfied  if  we  choose  M  =  k2 
Using  a  similar  argument  as  in  deriving  (3.16),  we  can  show  that  the  method 
represented  by  (3.19)  is  Ag-stable  if 


or 
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2  k-3 

D  1  2K  J  +  1  Q.E.D. 

Remark 

From  (3.18),  the  magnitude  of  the  error  constant  of  the  method 
2  k-3 

represented  by  (3.19)  when  D  =2+1  is  of  order 


2 

0((k-l )2k  /2) 


Corollary  3.8.2 

For  any  ae(0,ir/2)  and  any  C  >  0,  there  exist  A(a)-stable  3-step 
methods  of  order  3  with  error  constant  -C. 

Proof 

Assume  that  ae(0,Tr/2)  and  C  >  0  are  fixed.  Consider  the  class  of 
third  order  methods  represented  by 


s(z)  =  z(z2  +  24Cz  +  R2) 

where  R  .>  max{l  ,12C}.  The  magnitude  of  the  error  constant  of  the  method  is 
C.  It  is  obvious  that 

r(z)  =  2 [u ( z )  +  -] 

3 


where,  from  (3.11), 


u (z )  =  z2  +  24Cz  +  R2 

Hence  the  root  condition  is  satisfied.  The  locus  r(iw)/s(iw)  can  then  be 
expressed  as 

r(iw)  2  1 

- =  — [1  +  - '] 

s(iw)  iw  3u(iw) 

We  first  find  a  uniform  upper  bound  for  1 / (3u ( i w) )  over  -»  <  w  = 
Ry  <  oo.  if  R  >  12v^C,  then 
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24C  —  12C  — 

lu(iRy)  I  =  R2[(l-y2)2  +  <— )2y2]2  i  24RC[1  -  (— )2]2  1 12^RC 

R  R 

Thus 

11 
I - l< - — 

3u ( i  Ry )  R(36v^C) 

uniformly  in  ye(-~,“).  Using  a  similar  argument  as  that  in  (3.16),  the 
method  is  A( a) -stable  if  R  satisfies 

1 

- <_  si n (-  -  ct) 

R(36v#C)  2 

in  addition  to  the  assumption  that  R  >.  max{l  ,12/2C}.  Q.E.D. 

Remark 

2 

Note  that  the  roots  of  r(z)  have  magnitude  R  +1/3.  This  implies 
that  as  a  ->  tt/2 ,  the  polynomial  p(?)/(c-l)  has  roots  close  to  1.  The 
corollary  simply  demonstrates  the  inadequacy  of  the  error  constant  as  a 
measure  of  accuracy  for  k  =  3. 

3.6  CONCLUDING  REMARKS 

As  remarked  in  section  1.4,  an  ODE  solver  could  have  for  each  or¬ 
der  k  a  family  of  formulas  parameterized  by  a.  We  have  to  know  what  a  we 
want  before  we  can  select  the  right  formula.  While  the  order  k  is  chosen 
automatically,  the  angle  a  would  be  either  supplied  by  the  user  or  also 
determined  by  the  code.  As  an  illustration,  for  specified  values  of  k  and 
a,  we  could  choose  among  the  class  of  methods  of  the  form 

s (z )  =  (z  +  d)(z  +  D)k_1 

the  method  whose  parameters  d  and  D  satisfy  (3.16)  (the  values  for  M  and  K 
can  be  found  using  a  similar  argument  as  in  the  proof  of  corollary  3.8.1), 
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The  method  will  then  be  of  order  k  and  A(a)-stable.  We  note  that  the  above 
class  of  methods  can  be  parameterized  by  the  two  parameters  d  and  D,  or  by 
the  parameter  D  alone,  in  which  case  d  can  be  a  fixed  constant  or  a  func¬ 
tion  of  D  (d  has  to  be  of  order  o{D)  so  that  (3.16)  can  be  satisfied  when  a 
-*  tt/2 ) .  However,  the  selection  of  the  formulas  for  various  ct  might  best  be 
done  through  some  table  residing  in  the  code.  The  simplest  such  ar¬ 
rangement  would  limit  the  order  depending  on  a. 
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CHAPTER  IV  MINIMAX 

4.1  INTRODUCTION 

In  the  next  chapter,  we  shall  be  concerned  with  investigating  how 
large  the  angle  of  absolute  stability  can  be  among  methods  having  error 
constants  of  the  same  magnitude.  Since  an  analytical  solution  is  difficult 
to  obtain,  we  shall  solve  it  using  numerical  methods.  If  we  restrict  our 
attention  to  A  -stable  methods  only,  then  from  (3.2),  a  method  (r,s)  is 

oo 

A(a)-stable  if 

Re  r(iw)/s(iw) 

sup  -  -  £  cot  a 

-oo<w<oo  | Im  r(iw)/s(iw) | 
r(iw)fO 

Finding  the  limitation  on  the  angle  a  among  Aro-stable  methods  having  error 
constants  of  the  same  magnitude  is  then  equivalent  to  solving  for  the  in- 
fimum  of  the  quantity  on  the  left  hand  side  of  the  above  inequality  over 
all  values  of  the  parameters  bQ,b1 . bk_-|  subject  to  appropriate  con¬ 

straints.  In  this  chapter,  we  discuss  a  numerical  search  procedure  for 
solving  constrained  minimax  problems.  The  formulation  of  the  above  problem 
into  a  constrained  minimax  problem  is  given  in  more  detail  in  Chapter  V. 

A  formulation  of  the  constrained  minimax  problem  is  stated  in 
section  4.2.  A  feasible  descent  algorithm  is  discussed  in  section  4.3,  the 
convergence  of  which  is  proved  in  the  subsequent  section.  Some  modifica¬ 
tions  of  the  algorithm  are  suggested  in  section  4.5. 

4.2  STATEMENT  OF  THE  PROBLEM 

Consider  a  continuous  real-valued  function  f(x,y)  defined  on  the 
domain  X'xY,  where  Y  is  a  nonempty  compact  set  in  Rm,  and  X'  is  an  open  set 
in  In  containing  as  a  subset  the  set  X  described  by  the  following  nonlinear 
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inequality: 

X  =  {x£JRn  |  g(x,u)  £  0  for  all  ueU}  (4.1) 

in  which  U  is  a  nonempty  compact  subset  of  Ip  and  g  is  a  continuous  real- 

valued  function  defined  on  RnxU.  The  problem  is  to  minimize  the  objective 

function 

*(x)  =  max  f(x,y)  (4.2) 

yeY 

over  the  closed  (not  necessarily  bounded)  set  X.  The  set  X  is  called  the 

feasible  region  or  constraint  set  and  any  xeX  is  a  feasible  point. 

The  approach  that  is  used  in  solving  most  constrained  optimiza¬ 
tion  problem  goes  as  follows:  first,  Lagrange  multiplier  theorems  descri¬ 
bing  the  necessary  conditions  that  an  optimal  solution  must  satisfy  are 
formulated,  then  an  iterative  search  is  devised  to  locate  points  which  sa¬ 
tisfy  the  necessary  conditions.  Those  points  are  usually  called  stationary 
points.  To  this  end,  we  require  that  the  functions  f(x,y)  and  g(x,u)  be 
continuously  Fr^chet  differentiable  (Luenberger,  1969,  p. 1 72 )  with  respect 
to  x  on  X'xY  and  RnxU,  respectively.  It  can  then  be  proved  that  the  func¬ 
tion  <(>  is  Gateaux  differentiable  (Luenberger,  1969,  p .  1 71 )  at  each  xeX', 
the  Gateaux  differential  of  <j>  at  x  in  the  direction  d  being  given  by 
(Dem'yanov  and  Malozemov,  1974,  p. 1 88) 

9<j>  3f  ■ 

— (x)  =  max  (— (x,y),d)  (4.3) 

ad  ysR(x)  ax 

where  (  ,  )  is  the  scalar  product  of  two  vectors  and 

R(x)  =  {^Y  |  <|>(x)  =  f (x,y)}  (4.4) 

A  stronger  statement  is  in  p. 1 91  of  the  same  reference: 
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3f 

<f)(x+hd)  =  <j>(x)  +  h  max  ( — (x,y),d)  +  o(h;d)  (4.5) 

yeR(x)  3X 

where  o(h;d)/h  0  as  h  -»•  0  uniformly  in  d,  |  |d|  |  =  1  ( |  j  1 1  is  the  Eu¬ 
clidean  norm). 

4.3  A  FEASIBLE  DESCENT  ALGORITHM 

As  pointed  out  in  the  previous  section,  we  first  state  a  theorem 
describing  stationary  points  of  $  on  X  (cf.  (4.1)  and  (4.2))  and  then  pro¬ 
pose  a  feasible  descent  algorithm  to  find  a  stationary  point.  The  fol¬ 
lowing  notation  will  be  needed,  in  which  xeln,  yeRm,  ugR*3: 

Q(x)  =  {ueU  |  g(x,u)  =  0} 
af 

H(x)  =  {— (x,y)  |  yeR(x)} 
ax 

39 

H'(x)  =  { — (x.u)  |  ueQ(x)} 
ax 

L(x)  =  conv  [H(x)  u  H' (x)] 

s  s 

conv  S  =  {  Z  a • Z .  I  Z.gS,  a . > 0 ,  1< i< S ,  E  a.=l) 
i=l  11  1  1_  --  i=l  1 

r 

cone  S  =  {  E  X.Z,  |  Z.es,  x.>0,  lii<r} 
i=l  1  1  1  1 

The  function  f(x,y)  is  said  to  be  binding  at  x  if  yeR(x).  Similarly,  the 
constraint  g(x,u)  is  said  to  be  active  at  x  if  ueQ(x). 

Dem'yanov  and  Malozemov  (1974,  pp. 191-196, 203-204)  derive  the 
following  necessary  conditions  for  the  optimal  solution: 

Theorem  4.1 

If,  in  addition  to  the  assumptions  stated  in  section  4.2,  we  as- 
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sume  that  for  each  ueU,  g(x,u)  is  convex  in  x,  and  that  the  Slater  condi¬ 
tion  is  satisfied,  i.e.,  there  exists  xe»n  such  that 

max  g(x,u)  <0 
ueU 

then  a  necessary  condition  for  the  function  <j>(x)  to  achieve  its  minimum  on 
X  at  . a  point  x*eX  is  that 

conv  H(x*)  n  cone  [-  H'(x*)]  4  0 

An  equivalent  necessary  condition,  but  computationally  easier  to 
test,  can  be  derived  using  the  same  idea  as  in  Dem'yanov  and  Malozemov 
(1974,  pp. 146-148  in  which  U  is  assumed  to  be  a  finite  set  only): 

Corollary  4.1 

The  necessary  condition  in  theorem  4.1  is  equivalent  to 

0eL(x*)  (4.6) 

Definition 

A  point  x*eX  which  satisfies  (4.6)  is  called  a  stationary  point 

of  on  X. 


We  proceed  to  describe  a  feasible  descent  algorithm  for  finding 
stationary  points  of  <j>  on  X.  For  any  feasible  point  x,.a  nonzero  vector 
del  is  said  to  be  a  feasible  direction  on  X  at  x  if  there  exists  a  scalar 
TT  >  0  such  that  x+hdeX  for  all  he[0,h]  (Zoutendijk,  1966).  A  necessary 
condition  for  a  nonzero  vector  d  to  be  a  feasible  direction  at  xeX  can  be 
shown  to  be 

3g 

( — (x,u),d)  <.0  for  all  ueQ(x)  (4.7) 

8X 

We  note  that  some  authors  call  d  a  feasible  direction  if  (4.7)  holds.  It 
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can  be  proved  that,  when  the  assumptions  of  theorem  4.1  are  satisfied,  the 
closure  of  the  set  of  feasible  directions  at  x  is  equal  to  the  set  of  deRn 
satisfying  (4.7)  (Dem'yan°v  and  Malozemov,  1974,  p. 1 99) .  We  call  a  nonzero 
vector  d  a  feasible  descent  direction  for  <|>  on  X  at  x  if  there  exists  a  po¬ 
sitive  scalar  h*  such  that,  for  any  he(Q,h*],  x+hdeX  and  (x+hd )  <  <f>(x) 
(Zoutendijk,  1966). 

The  feasible  descent  algorithm  works  as  follows.  A  feasible 
starting  point  Xg  is  given.  Suppose  that  we  have  obtained  points 
x0,x-,,...,xveX.  If  xv  is  a  stationary  point,  the  algorithm  terminates. 
Otherwise,  a  new  point  xv+1eX  is  determined  such  that  4>(xv+1)  <  <j>(xv),  or 
it  is  concluded  that  no  stationary  point  exists.  The  latter  may  occur  sin¬ 
ce  X  is  not  bounded.  The  determination  of  xv+1  consists  of  two  parts:  a 
feasible  descent  direction  dv  at  xy  is  found,  then  a  step  is  made  in  the 
direction  dv  by  choosing  a  step  size  hv  so  that  the  new  point  xv+1  = 
xv+hvdv  is  in  X  and  satisfies  <f>(xv+])  <  <fr(xv).  If  the  tw0  conditions  can 
be  satisfied  by  any  h  >  0,  or  if  ||x  ||  becomes  unbounded  as  v  increases, 
then  we  conclude  that  no  stationary  point  exists. 

The  feasible  descent  directions  can  be  determined  using  the  bin¬ 
ding  functions  and  the  active  constraints.  Since  the  function  <j>(x)  is  not 
Frdchet  differentiable,  such  strategy  could  result  in  convergence  to  a  non¬ 
stationary  point  (example  in  Dem'yanov  and  Malozemov,  1974,  pp. 76-82  where 
X  =  in),  Such  a  phenomenon  is  known  as  zigzagging.  It  occurs  when  the 
sequence  of  steps  becomes  progressively  very  small  not  because  a  stationary 
point  is  in  the  vicinity  but  because  new  binding  functions  or  new  active 
constraints  are  encountered.  One  remedy  is  to  take  into  consideration  not 
just  the  binding  functions  and  the  active  constraints  but  also  the  "nearly" 
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binding  functions  and  the  "nearly"  active  constraints.  This  also  serves  to 
prevent  binding  functions  and  active  constraints  from  "disappearing"  due  to 
roundoff  errors. 

Let  e,y  be  nonnegative  numbers  and  define  the  sets 

Re(x)  =  {yeY  |  <f>(x)  -  f(x,y)  <  e}  (4.8) 

Qy(x)  =  {ueU  |  -  y  <  g(x,u)  <0} 

3f 

H  (x)  =  { — (x,y)  |  yeR  (x)} 

3X  £ 

3g 

H'(x)  =  {— (x,u)  |  ueQ  (x) > 

M  3X  U 

Ley(x)  =  conv  [He(x)  i  u  H^(x)]  (4.9) 

Note  that  RQ(x)  =  R(x),  Qq(x)  =  Q(x),  HQ(x)  =  H(x),  HJ,(x)  =  H'(x),  and 
Lqq(x)  =  L(x).  Moreover,  the  sets  H£(x)  and  H^(x)  are  compact  and  hence 

Lr(x)  is  a  compact  convex  set.  Let  z  (x)  be  the  point  of  L  (x)  nearest 

ey  ey 

the  origin,  i  .e. , 


|z  (x) 

1  ey '  ' 


mm 

zgL  (x) 

ey 


If  z  (x  )  =  0,  or  equivalently,  OeL  (x  ),  then  we  call  x  an  (e,y)- 

t-  (J  v  ^  l-l  V  V 

quasistationary  point  of  <j>  on  X,  in  which  case  e  and  y  are  reduced  and  the 
search  continues  with  the  new  e  and  y.  Note  that  (0,0)-quasistationarity 
is  equivalent  to  stationarity.  If  z£y(xv)  f  0,  define 


d  =  - 

V 


z  (x  ) 

ey'  v' 

Z  (X  ) 

ey'  v 


(4.10) 


The  vector  d  is  a  feasible  descent  direction  as  shown  by  the  following 


1 emma : 
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Lemma  4.1 

For  any  e  >  0  and  y  >  0,  if  xv  is  not  (e.y)-quasistationary,  then 
dv  as  defined  in  (4.10)  is  a  feasible  descent  direction  for  <f>  on  X  at  xv. 
Furthermore, 

9<j> 

—  (x  )  <  -  |  |z  (x  )  I  I 
j  '  v7  —  ey  v  1  1 

v 

Proof 

Since  z  (x  )  is  the  point  of  L_„(xM)  nearest  the  origin,  it  can 

ey  v  ey  v 

be  shown  that  (Dem'yanov  and  Malozemov,  1974,  p.252) 

IUEp(xp)l|2£(z,z£p(xp))  for  all  zeLEp(xp) 

Thus,  from  (4.10),  we  have 

(z,dv)<-  ||zey(xv)||<0  for  all  zeLE]j(xv)  (4.11) 

In  particular, 

(— (xv,y),dv)  i-  ||zey(xv)||  for  all  yeR(xv) 

3x 

Hxv,u),dv)  1-  llzey(xv)M  for  all  ueQ(xv) 

3x 

Using  (4.3)  and  (4.5),  and  a  similar  argument  for  the  function 

max  g(x,u) 

U€U 

it  is  then  obvious  that  dv  is  a  feasible  descent  direction.  Q.E.D. 

After  determining  the  feasible  descent  direction  dv,  we  next  want 

to  take  a  step  in  that  direction.  The  stepsize  hv  can  be  found  by  the  fol¬ 
lowing  step  selection  rule,  which  is  a  modification  of  the  Armijo  rule  (Po- 
lak,  1971,  p.36): 
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Let  y  >  0,  ge(0tl),  and  ae(0,l)  be  fixed  constants.  The  values  y 
=  1,  ge(0.5,0.8),  and  a  =  0.5  are  recommended  by  Polak  (1971,  p.36).  We 
start  with  an  initial  value  for  hv  given  by 


If  hv  leads  to  an  infeasible  point,  it  is  reduced  by  the  factor  3  until  the 
feasibility  constraint  is  satisfied.  Note  that  since  d  is  a  feasible 

V 

direction,  at  most  a  finite  number  of  reductions  is  needed.  There  are  then 
two  possibilities.  If 

1  4>(xv)  -  °hj  |zep(xv)|  |  (4.12) 


then  we  define  x  , 7  =  x  +  h  d  .  Otherwise,  the  stepsize  h  is  reduced 
successively  by  the  factor  3  until  the  above  inequality  is  satisfied.  The 
hv  corresponding  to  the  smallest  <j>(xv+hvdv)  obtained  in  this  iteration  is 
the  stepsize  used  in  defining  xv+1 ,  i.e.,  xv+1  =  xv  +  hydv.  Since  d^  is  a 
feasible  descent  direction,  the  number  of  reductions  must  be  finite.  In 
case  llxv+hx)dv||  is  large,  it  is  conceivable  that  no  stationary  point  ex¬ 
ists,  hence  the  algorithm  terminates. 

The  constant  y  ensures  that  the  initial  stepsize  is  not  too 
small,  whereas  the  choice  hv_-|/3  enables  us  to  try  a  possible  larger  step. 
We  remark  that  the  original  Armijo  rule  uses  y  as  the  initial  trial  and 
chooses  as  the  stepsize  the  first  h  such  that  x  +h  d  eX  and 

V  V  V  V 


3<J) 

<f>(x  +h  d  )  <  <j>(x  )  +  ah - (x  ) 

'  v  v  v'  —  T '  v  v*  .  v 

3d 

v 


Since  <f>(x)  is  not  Frechet  differentiable,  we  replace  3<|>(xv)/3d  by 
-  I lzey(xv)| |  which  is,  from  lemma  4.1,  greater  than  the  former  quantity. 
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Note  that  we  can  choose  as  the  stepsize  the  first  hv  such  that  xv+hvdveX 
and  (4.12)  is  satisfied.  The  convergence  theorem  proved  in  the  next  sec¬ 
tion  still  holds,  and  in  fact  holds  for  any  line  minimization  rule  which 
gives  a  better  point,  i.e.,  smaller  4>(x  +i ) »  than  the  rule  just  mentioned. 
Since  it  is  more  desirable  to  have  a  larger  decrease  along  the  descent 
direction,  we  choose  as  the  stepsize  the  one  corresponding  to  the  smallest 
d>(x  +h  d  )  among  the  values  tried. 

Suppose  that  stationary  points  exist.  For  a  given  pair  of  num¬ 
bers  (e,y),  the  direction  and  step  selection  algorithm  as  described  could 
be  repeated  until  an  (e.y)-quasistationary  point  is  found.  In  practice,  we 
need  only  continue  the  iteration  until  a  point  xv^X  is  found  such  that 
|  |zey(xv) | |  <  p  for  some  positive  number  p.  Then  the  three  numbers  e,y,p 
are  reduced  and  the  algorithm  is  applied  again  using  the  new  e,y,p,  and  the 
current  xv  as  an  initial  point.  The  entire  algorithm  will  be  referred  to 
as  the  M-algorithm  and  the  portion  in  which  e,y,P  are  fixed  will  be  refer¬ 
red  to  as  the  (e,y,p)-algorithm.  Let  {ek  |  k  >.  1  >,  {yk  I  k  >.  1  >,  and  Tp k I 
k  >.1}  be  the  strictly  decreasing  sequences  of  positive  numbers  used  by  the 
M-algorithm.  We  use  the  superscript  k,  i.e.,  x  ,  to  denote  the  point  re¬ 
turned  by  the  (e,y,p)-algorithm  with  e  =  ek,  y  =  uk,  p  =  p^.  In  other 

k  .  . 

words,  each  x  ,  k  >_  1 ,  satisfies 

I  lzvk(xk)  I  I  <  pk 

The  points  generated  by  the  (e,y,p) -algorithm  for  the  values  ek,yk,Pk  is 

k  k  k-1  i 

denoted  by  the  additional  subscript  v,  i.e.,  x\  Note  that  xQ  =  x  .  The 

k  k 

fact  that  there  exists  a  finite  v.  such  that  x  =  x„  can  be  shown  using 

K  vk 

theorem  4.2  (cf  section  4.4).  The  convergence  of  the  sequence  {xk  |  k  >_ 1 } 
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to  a  stationary  point  is  discussed  in  the  next  section. 

We  finally  remark  that  an  approximation  for  the  vector  z  (x)  can 

be  obtained  as  follows.  Discretize  the  sets  R_(x)  and  Q,,(x)  and  define  the 

c  y 

matrix  D  whose  columns  consist  of  vectors  in  Hp(x)  and  H'(x)  evaluated  at 

^  y 

the  discretized  points  (cf.  (4.9)).  Then  z£y(x)  can  be  approximated  by  Da* 
where  a*  is  the  optimal  solution  of  the  following  convex  quadratic  program¬ 
ming  problem: 

minimize  aWDa 

subject  to  a.  >0,  1»l,...,q  a  =  [a-, ,«2, . . . ,aq]T 
a-|+a2+. .  *+aq  =  1 

The  above  problem  can  be  solved  by  a  number  of  methods,  e.g.,  the  "first 
symmetric  variant  of  the  simplex  method"  (Van  De  Panne,  1975,  pp. 270-280) 
which  is  shown  to  converge  to  the  optimal  solution  in  a  finite  number  of 
iterations. 

4.4  CONVERGENCE 

In  this  section,  we  prove  that  the  M-algorithm  described  in  the 
previous  section  converges  to  a  stationary  point  if  it  exists.  The  proof 
is  in  two  stages:  the  convergence  of  the  (e,y,p)-algorithm  for  positive 
e,y,p  and  then  the  convergence  of  the  M-algorithm  to  a  stationary  point  as 
ek  *  yk  an<^  pk  0*  The  existence  of  stationary  points  is  assured 

k  i 

if  the  sequence  (xv  I  0  <  v  <  vk,  k  >  1 }  generated  by  the  M-algorithm  is 
bounded  or  the  set 

X(x°)  =  {x^X  |  4>(x)  l<Kx°)) 

is  bounded,  where  x°gX  is  a  given  starting  point  for  the  M-algorithm.  For 
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simplicity,  we  assume  that  the  latter  holds. 

Lemma  4.2 

For  any  e  >  0,  there  exists  6  =  6(e)  >  0  such  that 

R{x')cR£(x) 

for  all  | | x 1 -x | |  <  6,  x',x£X(x^). 

Proof 

The  assertion  follows  directly  from  the  definition  of  R(x)  and 
Re(x)  (cf.  (4.4)  and  (4.8))  and  from  the  uniform  continuity  of  <f>(x)  on 
X(x°)  and  f(x,y)  on  X(x°)xY.  Q.E.D. 

Lemma  4.3 

For  any  e  >  0,  n  >  0,  there  exists  h^  =  h^U.n)  >  0  such  that 

4>  (x+hd)  -  4>  ( x )  3f 

- max  ( — (x,y),d)  <  0 

h  yeRe(x)  9x 

for  all  he(0,h(j)],  xeX(x°),  and  d€Kn  such  that  x+hdeX(xD)  and  ||d||  =  1. 
Proof 

From  the  uniform  continuity  of  9f/9x  on  X(x°)xY,  there  exists  po¬ 
sitive  h'  =  h'(n)  such  that 

9f 

f  (x+hd,y)  =  f(x,y)  +  h(— (x,y),d)  +  hA(x,y,hd) 

9X 

where  |A(x,y,hd)|  <  n  for  any  xgX(x^),  y^Y,  he(0,h'],  and  deKn  such  that 
x+hdeX(x^)  and  | |d| |  =  1 . 

By  lemma  4.2,  there  exists  positive  h"  =  h"(e)  such  that 

R(x+hd)  cRjx) 

for  any  xgX(x°),  he(0,h"],  and  de$n  such  that  x+hdeX(x°)  and  ||d||  =  1. 
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Let  h^  =  min{h',h"}.  Then  for  any  xeX(x^),  he(0,h^],  and  detn 

such  that  x+hdeX{x^)  and  ||d||  =  1, 

4)  (x+hd )  =  max  f(x+hd,y)  _<  max  f(x+hd,y) 
yeR(x+hd)  y6R£(x) 

3f 

<  max  f(x,y)  +  h  max  ( — (x,y),d)  +  hn 
yeR£(x)  yGR£(x)  3x 


3  f 

=  <|>(x)  +  h  max  ( — (x,y),d)  +  hn  Q.E.D. 

y^£(x)  3  x 


Theorem  4.2 

For  any  e  >  0  and  y  >  0,  if  (xv  |  v  0}  is  the  sequence  of 
points  generated  by  the  (e,y,0)-algorithm  with  initial  value  x0gX(x°),  then 

I  I  Zey  (Xv )  I  I  0  as  V  +  »,  v  >_  0 

Proof 

Let  n  >  0  be  given.  Define  (cf.  lemma  4.3,  4.4) 
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h*  =  minth^e.il-aJnJ.h^Cy.rvJ.y} 

Since  (4>(xv)  |  v  >  0}  is  a  monotone  sequence  bounded  below  by 

inf  <Hx)  =  1nfn  <f>(x)  >  -  00 
xeX  xgX(x  ) 

there  exists  a  positive  integer  N  =  N(e,y,n)  such  that 

<t>(xv)  -  4>(xv+i )  <  a3h*n  for  all  v  >  N 

Let  hv  be  the  stepsize  last  tried  by  the  step  selection  rule  (not  neces¬ 
sarily  the  stepsize  chosen).  Then  xv+hydyeX(x^)  and 

<f>(xv)  -  <Kxv+1)  >  +(xv)  -  4>(xv+hvdv)  >  ahv||zey(xv)|| 

Consider  v  11.  If  hy  3h*,  then 

<f>(xj  -  *(x  +1)  agh*n 
1 1  z  (x  )  < — - - —  < - <n 

ey '  v '  ”  ^  ^ 

Suppose  that  h  /3  <  h*.  Since  h  <  y,  there  are  two  possibilities: 


Case  1. 


♦(ypjyttg  -  »og  > 

hv/6  "  v 


Since,  from  (4.11), 


(— (xv.y),d  )  1  - 

9X 

using  lemma  4.3,  we  obtain 


ZE„(XJ 


for  all  yeR£(xv) 


z  (x  ) | |  <  n 
1  ey '  v  1  1 


Case  2. 


^(xv+[Rv/g]dv)  >  0  >  *(xv) 
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or 

<l<(x  +[h  /g]d  )  -  *(x  ) 

V  V  V  ^  >  0 

fiv/8 

Using  (4.11)  and  lemma  4.4,  we  have 

IIWVH  <  n  Q.E.D, 

One  immediate  corollary  of  the  above  theorem  is  that  the 
(e,v,p )-algorithm  terminates  after  a  finite  number  of  iterations  if  e,y,p 
are  positive.  Note  that  the  above  theorem  holds  if  the  new  point  is  de¬ 
fined  by  fi^,  i.e.,  +  h^d^,  as  remarked  in  the  previous  section. 

To  prove  that  the  M-algorithm  yields  stationary  points  of  <j>  on  X, 
we  have  to  show  that  for  small  enough  ek»yk»pk»  if  x  is  close  to  an 
(£k'yk^"quasistationary  P°int>  then  xk  also  dose  to  a  stationary  point. 
We  shall  need  the  following  lemmata: 

Lemma  4.5 

For  any  and  any  n  >  0,  there  exist  positive  6  =  <s(n,x),  7  = 
e(n,x),  and  y  =  y(n,x)  such  that  the  following  hold  for  any  ||x'-x||  <  S, 
x'ex: 

(i)  For  any  ee[0,7],  every  element  of  R£(x')  is  within  an  n-neighborhood  of 
an  element  of  R(x). 

(ii)  For  any  ye[0,y],  every  element  of  Qy(x')  is  within  an  n-neighborhood 
of  an  element  of  Q(x). 

Proof 

Suppose  there  exist  a  positive  number  n',  a  positive  sequence 
|  i  >  01  converging  to  zero,  and  a  sequence  of  points  (x^eX  |  i  >_ 0>  con- 


verging  to  x  such  that  for  any  i  >_  0,  there  exists  y^R^  (x..)  which  is  not 

in  an  n '-neighborhood  of  any  \^R(x),  i.e., 

Il^i  -  v|  |  >_  r\ 1  for  all  v€R(x)  (4.13) 

Since  y..eY  for  all  i  >  0,  there  is  a  subsequence  {y^  |  iel},  I  c  {0,1,...}, 
converging  to  y£Y.  For  each  iel, 

i>(xi)  -  f (xi  ,yi )  1  e.j 

By  continuity,  as  i  ->  ®,  iel, 

4>(x)  -  f(x,y)  <0 

which  implies  that  yeR(x).  It  follows  that  for  sufficiently  large  i,  there 
exists  yeR(x)  such  that  ||yi-y||  <  n '  contradicting  (4.13). 

To  prove  (11'),  we  have  to  distinguish  between  two  cases:  when 
Q(x)  =  0  and  when  Q(x)  +  0.  The  first  case  is  trivial  while  the  latter  can 
be  proved  using  a  similar  argument  as  above.  Q.E.D. 

Lemma  4.6 

Let  G  be  a  compact  set  in  Kn  and  zQ  be  the  point  of  conv  G 
nearest  the  origin.  Suppose  that  ||zg||  0.  If  G1  is  a  compact  set  such 

that  every  element  of  G'  is  within  an  1 1 zQ |  |/(2v/n)-neighborhood  of  an 
element  of  G,  then  for  any  z'econv  G', 

Hz' 1 1  >  l  |z0ll/2 

Proof 

From  the  hypotheses,  it  can  be  proved  that  every  element  of 
conv  G'  is  within  an  | |zQ| | /2-neighborhood  of  an  element  of  conv  G,  i.e., 
for  any  z'econv  G‘,  there  exists  zeconv  G  such  that 
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It  follows  immediately  that 

ll2'll  >  llz0M  -  llz0ll/2  =  llz0li/2  Q.E.D. 

Lemma  4.7 

For  any  xeX,  there  exist  positive  6  =  6(x),  7  =  7(x),  and  7  = 

7(x)  such  that  for  all  ||x'-x||  <  8,  x'eX,  ee[0,7],  ue[0,7], 

l|zey(x')ll  >  ■||z00(x)||/2 

Proof 

The  proof  for  the  case  when  ||z00(x)||  =  0  is  trivial. 

Suppose  that  ||zqq(x)||  ^  0,  i.e.,  x  is  not  a  stationary  point. 
From  the  uniform  continuity  of  the  functions  9f/9x  and  9g/9x  in  yey  and 
u€U,  respectively,  we  can  find  n  =  n(x)  >0  such  that  the  following  hold 
for  any  |  |x'-x| |  <  n,  x'eX: 

9f  I lznn(x) | | 

1 1 — (x,y) - (x1  ,y‘ )  1 1  <  — w3— —  for  all  y.y'eY,  |  |y-y'  1 1  <  n 

9x  9x  2/n 

99,  ,  ||znn(x)|| 

1 1 — (x,u)  -  — (x1  ,u')|  |  <  - — -  for  all  u,u'ell,  I  I u-u '  I  I  <  n 

ax  3x  2v4T 

Applying  lemma  4.5,  we  obtain  positive  numbers  s'  =  S'(x),  7  =  7(x),  and  7 

=  y(x)  such  that  the  hypotheses  of  lemma  4.6  are  satisfied  with 

9f  9g 

G  =  ( — (x,y)  |  yeR(x)}  u  {— (x,u)  |  ueQ(x)) 

9X  9X 

9f  9g 

G'  =  { — (x',y')  |  y'eR  (x')>  u  {— (x',u')  |  u'eQ  (x')> 

9X  9X  y 

for  any  x'eX,  ||x'-x||  <  6  =  min{n,<5'},  ee[0,7],  and  ye[0,7].  It  follows 

from  lemma  4.6  that 


l|zey(x')ll  1  I !z00(x) 1 1/2 


Q.E.D. 
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Theorem  4.3 

Every  limit  point  of  the  sequence  (xk  |  k  1}  generated  by  the 
M-algorithm  is  a  stationary  point  of  $  on  X. 

Proof 

Let  x*  be  a  limit  point  of  {xk  |  k  >_  1}.  Then  there  is  a  subse¬ 
quence  (xk  |  kel},  I  c  {0,1,...},  converging  to  x*.  Using  lemma  4.7,  there 
exists  a  positive  integer  K  such  that  for  all  k  >_  K,  kel, 

I  |Zqq(x*)  II/2  1  I  l^ky^  ^  ^  pk 

Since  pk  -*  0  as  k  ->  °°,  so  must 

1 1 z0Q(x*)  I  i  0  Q.E.D, 

Theorem  4.3  does  not  guarantee  the  convergence  of  the  sequence 
(xk  |  k  >  1>.  If  we  further  assume  that  the  function  <j>  has  at  most  a  finite 
number  of  stationary  points,  it  could  be  shown  that  the  generated  sequence 
will  converge  to  a  stationary  point  of  <J>  on  X. 

4.5  MODIFICATIONS 

Since  X  is  described  by  a  nonlinear  inequality,  the  function 
g(x,u)  can  be  defined  arbitrarily  in  the  sense  that  the  same  region  is  ob¬ 
tained  if  we  replace  g(x,u)  by  g(x,u)  times  a  positive  function.  Subse¬ 
quently  the  choice  of  feasible  descent  directions  is  also  arbitrary  depen¬ 
ding  on  how  g(x,u)  is  defined.  Such  an  arbitrariness  can  also  be  seen  from 
theorem  4.1.  The  vectors  in  H'(x*)  can  be  adjusted  by  any  positive  number 
and  the  necessary  condition  for  x*  to  be  a  global  minimum  still  holds  with 
H'(x*)  replaced  by 
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39 

H'(x*;t)  =  {t(x*,u) — (x*,u)  |  ueQ(x*)} 

3x 

where  t(x*,u)  >  0  for  all  ueQ(x*).  The  convergence  property  of  the  M- 
algorithm  will  not  be  altered  if  we  require  that  the  numbers  t(x,u)  be  po¬ 
sitive  and  be  bounded  above  by  a  constant.  One  particular  choice  for 
t(x,u)  is  motivated  by  the  following  consideration  in  the  two  dimensional 
space: 

Suppose  that  in  choosing  a  feasible  descent  direction  we  approxi¬ 
mate  the  constraint  set  by  a  disc  and  use  a  linear  approximation  for  the 
objective  function.  The  resulting  problem  is  to  minimize 

f(x)  =  cTx  c.xgk2 
over  the  constraint  set  defined  by 

X  =  (xeK2  |  g (x)  =  (x  -  a)T(x  -  a)  -  r2  £  0),  aeR2 
Note  that  with  appropriate  scaling  of  the  decision  variables,  a  convex  con¬ 
straint  set  can  be  approximated  by  a  disc.  It  will  be  seen  that  the  radius 
of  the  disc  is  immaterial  in  the  determination  of  feasible  descent  direc¬ 
tion. 

Suppose  that  we  start  with  an  initial  point  x°  lying  on  the  boun¬ 
dary  of  X  (see  figure  4.1).  By  the  linearity  of  the  objective  function, 
the  optimal  solution  x*  must  lie  on  the  boundary.  Thus  the  vector  x*-x° 
gives  a  feasible  descent  direction  which  leads  to  the  optimal  point  in  one 
step.  Geometrically,  x*-x°  is  the  vector  which  bisects  the  angle  formed  by 
the  vectors  af(x°)/ax  and  3g(x°)/ax.  Using  the  Kuhn-Tucker  theorem  (Luen- 
berger,  1969,  p.249),  we  can  show  that  x*-x°  is  given  by 
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/ 


f(x)  -  f(x°) 

Figure  4.1.  A  Two  Dimensional  Case 


x*  -  x®  =  - 


r  3f  n  1 1 9f (xu)/ax  [  |  3g  0 

- C — — (x  )  +  - 0 - ; - (x  )> 

[cl  I  3x  |  |  3g(xU)/3X  I  I  3X 


|  [c  I  1  3x  |  |  3g(x  )/  3X  |  |  3X 

Note  that  the  direction  of  x*-x°  is  independent  of  r,  the  radius  of  the 
disc.  Thus  it  seems  like  it  may  be  advantageous  to  scale  3g(x)/3x  so  that 
it  has  the  same  norm  as  the  vector  3f(x)/3x.  Based  on  this  observation,  we 
could  pick  t(x,u)  so  that  the  vectors  t(x,u) 3g(x,u)/3x,  ueQ^(x),  have  norm 
equal  to  the  smallest  norm  of  the  vectors  3f(x,y)/3x,  y=R£(x). 

The  above  example  also  indicates  that  an  appropriate  scaling  of 
the  variable  x  is  desirable.  One  common  way  is  to  scale  x  by  a  linear 
transformation  matrix  T,  i.e., 


x  =  TC 


Using  linear  algebra,  we  can  show  that  the  above  scaling  is  equivalent  to 
the  algorithm  in  which  the  next  point  is  defined  by 
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x 


v+1 


x  +  h  TTTd 

V  V  V 


instead  of  x  , ,  =  x  +  h  d  .  We  shall  not  go  into  details. 

V+l  V  V  V  3 


4.6  CONCLUDING  REMARKS 

We  have  described  an  algorithm  for  finding  only  stationary  points 
of  *  on  X.  A  stationary  point  need  not  be  a  local  minimum,  it  can  also  be 
a  local  maximum  or  a  local  saddle  point.  The  first  possibility  can  only 
occur  at  the  starting  point  x°  if  there  exists  a  yeR(x°)  such  that 
3f(x°,y)/8x  =  0.  The  second  possibility  can  be  tested  by  exploring  dif¬ 
ferent  directions  to  see  if  <p  can  be  further  minimized.  Some  sufficient 
conditions  for  a  point  x*  to  be  a  local  minimum  are  suggested  by  Dem'yanov 
and  Malozemov  (1974,  pp. 125-126).  Note  that  if  the  convexity  or  the  Slater 
assumption  (cf.  theorem  4.1)  fails  to  hold,  the  M-algorithm  can  still  be 
used  to  locate  local  minima  of  <)>  on  X  provided  that  the  above  tests  are 
performed  on  the  generated  limit  points. 

If  the  function  <j>  is  not  convex,  there  may  be  more  than  one  local 
minimum,  and  in  the  best  case,  we  shall  only  find  one  of  them.  A  partial 
remedy  is  to  apply  the  M-algorithm  several  times  using  different  starting 
points. 

It  should  also  be  remarked  that  the  formalism  is  general  enough 
to  include  cases  when  X  is  defined  by  more  than  one  inequality  constraint. 
However,  because  of  the  Slater  condition,  problems  involving  equality  con¬ 
straints  can  only  be  handled  if  the  equality  constraints  can  be  solved  to 
reduce  the  original  problem  to  a  new  problem  of  smaller  dimension  having 
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only  inequality  constraints. 
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CHAPTER  V  ACCURACY  VS  A(a)-STABILITY 

5.1  DESCRIPTION  OF  THE  PROBLEM 

In  this  chapter,  we  are  concerned  with  discovering  the  upper 
bound  on  the  angle  of  absolute  stability  (definition  in  section  1.2)  for 
methods  in  L[k]  having  a  given  error  constant.  For  each  A  >  o,  define 
A*(a)  to  be  the  supremum  of  all  a  for  which  there  exists  an  A(a)-stable 

i/ 

method  having  error  constant  -A  ,  and  be  zero  if  no  A(0)-stable  method  with 

that  error  constant  exists.  Note  that  the  problem  independent  part  of  the 

leading  term  of  the  asymptotic  error  expansion  (2.1 )  of  such  methods  has 
k 

magnitude  ( h a)  .  Possible  values  for  A*(A)  are  shown  in  figure  5.1.  Our 


Figure  5.1.  Possible  values  for  A*(a) 

objective  is  to  find  A* (a)  for  A  >  0.  We  expect  that  in  all  probability 
A*(a)  is  strictly  increasing  on  some  semi-infinite  interval  in  (0,«). 

From  sections  3.3  and  3.5,  the  A*( a)  values  for  A  >  0  when  k  = 
1,2,3  are  as  shown  in  figures  5.2,  5.3,  and  5.4,  respectively. 

An  analytical  solution  for  the  general  case  when  k  ^4  is  diffi¬ 
cult.  We  thus  find  the  A*(a)  values  numerically.  A  mathematical  formula- 
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tion  is  discussed  in  the  next  section. 

We  remark  that  the  results  from  section  3.5  provide  a  "lower 
bound"  for  A*(a)  over  some  semi-infinite  interval  in  (0,~).  Moreover, 

A*(a)  approaches  the  line  a  =  ir/2  as  A  »  at  least  as  fast  as 
0(A‘2k/(k-2)).  It  also  follows  from  section  3.4  that  A*(a)  =  0  for  0  <  a  < 
l/[2(3k)'/k], 

5.2  MINIMAX  FORMULATION  OF  THE  PROBLEM 

A  mathematical  formulation  of  the  problem  described  in  the  pre¬ 
vious  section  is  suggested  by  condition  (3.2).  To  use  that,  we  consider 
only  methods  which  are  A  -stable.  For  reasons  that  will  be  apparent  later, 

oo 

we  further  restrict  our  attention  to  methods  that  satisfy  the  strict  root 
condition.  Since  convergent  methods  satisfy  the  root  condition  (cf.  sec¬ 
tion  1.3),  and  from  theorem  3.5,  A(0)-stable  methods  are  almost  A^-stable 
save  for  possible  simple  roots  on  the  imaginary  axis,  the  above  restric¬ 
tions  lead  to  only  a  "slightly  smaller"  feasible  region.  We  note  that 
there  exist  methods  which  are  not  stable  but  whose  locus  r(iw)/s(iw), 
v^(-OT,00),  is  outside  the  wedge  S&  for  some  a  >  0,  e.g.,  the  5-step  method 
represented  by 

(bg  ,b-|  ,b2  ,bg  ^ )  =  (0,15, .14, 3,. 007 ) 

From  (3.2),  an  /^-stable  method  satisfying  the  strict  root  condi¬ 
tion  is  A(a)-stable  if  and  only  if 
Re  r(iw)/s(iw) 

_  - cot  a  for  all  w _>  0  (5.1 ) 

|  Im  r(iw)/s(iw)j 

We  need  only  consider  the  semi-infinite  interval  we[0,“)  because  the  com¬ 
plex  conjugate  of  r(iw)/s(iw)  is  r(-iw)/s(-iw)  so  that  the  rational  func- 
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tion  is  discussed  in  the  next  section. 

We  remark  that  the  results  from  section  3.5  provide  a  "lower 
bound"  for  A*(a)  over  some  semi-infinite  interval  in  (0,e°).  Moreover, 

A*(a)  approaches  the  line  a  =  u/2  as  A  at  least  as  fast  as 
0(A-2k/(k“2)).  It  also  follows  from  section  3.4  that  A*(a)  =  0  for  0  <  a  < 
l/[2(3k),/k]. 

5.2  MINIMAX  FORMULATION  OF  THE  PROBLEM 

A  mathematical  formulation  of  the  problem  described  in  the  pre¬ 
vious  section  is  suggested  by  condition  (3.2).  To  use  that,  we  consider 
only  methods  which  are  A  -stable.  For  reasons  that  will  be  apparent  later, 

co 

we  further  restrict  our  attention  to  methods  that  satisfy  the  strict  root 
condition.  Since  convergent  methods  satisfy  the  root  condition  (cf.  sec¬ 
tion  1.3),  and  from  theorem  3.5,  A(0)-stable  methods  are  almost  Aro-stable 
save  for  possible  simple  roots  on  the  imaginary  axis,  the  above  restric¬ 
tions  lead  to  only  a  "slightly  smaller"  feasible  region.  We  note  that 
there  exist  methods  which  are  not  stable  but  whose  locus  r(iw)/s(iw), 
weC-00,00),  is  outside  the  wedge  Sa  for  some  a  >  0,  e.g.,  the  5-step  method 
represented  by 

(bQ,bi  ,b2,b3,b4)  =  (0,15, .14, 3, .007) 

From  (3.2),  an  /^-stable  method  satisfying  the  strict  root  condi¬ 
tion  is  A  (a) -stable  if  and  only  if 
Re  r(iw)/s(iw) 

-  - £  cot  a  for  all  w £  0  (5.1 ) 

|  Im  r(iw)/s(iw)| 

We  need  only  consider  the  semi-infinite  interval  we[0,°°)  because  the  com¬ 
plex  conjugate  of  r(iw)/s(iw)  is  r(-iw)/s(-iw)  so  that  the  rational  func- 
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tion  on  the  left  is  synmetric  about  w  =  0.  To  indicate  explicitly  the  fact 
that  both  r(z)  and  s(z)  are  actually  functions  of  b  =  (b^b^ .... *bjc_i )  and 
z,  we  write  them  as  r(b,z)  and  s{b,z),  respectively.  Since  |s(b,iw)|  >  0 
for  any  w  1  0,  (5.1)  is  then  equivalent  to 

Re  r(b,iw)s(b,-iw) 


sup  -  — - : — - 1  cot  a 

wlO  | Im  r(b,iw)s(b,-iw) j 


(5.2) 


The  restriction  that  the  method  satisfies  the  strict  root  condition  ensures 
that  the  numerator  and  denominator  of  the  rational  function  on  the  left  are 
not  zero  simultaneously;  so  does  the  Aro-stability  assumption. 

Suppose  that  a  value  of  A  (>  0)  is  given.  We  then  solve  for  the 
infimum  of  the  quantity  on  the  left  of  the  inequality  (5.2)  over  the  family 
of  methods  in  i[ k]  which  are  A^-stable,  satisfy  the  strict  root  condition, 
and  have  error  constant  of  magnitude  A^.  Hopefully,  the  arccotangent  of 

the  infimum  is  the  value  of  A*(a). 

We  proceed  to  put  the  problem  into  the  form  described  in  section 

4.2.  As  remarked  in  section  4.6,  the  equality  constraint 

1  k  b . 


nr  E  A 

2K  j=0  j+1 
j  even 

can  be  used  to  reduce  the  dimension  of  the  problem  by  one.  One  way  is  to 
use  b1,b2,...,bk_1  as  the  decision  variables  and  set 

b  .  ,k.l  kE  bj_ 

0  2k  j=2  j+1 

j  even 

We  denote  the  decision  variables  by  x  =  (xltx2,...  »x|<_-|)  where  Xj  =  b..,  j 
=  1,2,...  ,k-l ,  and  write  r(x,z),  s(x,z)  for  r(b,z),  s(b,z),  respectively, 
implying  that  the  equality  constraint  has  been  eliminated. 
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The  feasible  region  should  be  the  set  of  xeSk-1  at  which  the  po¬ 
lynomials  (in  z)  r(x,z)  and  s(x,z)  have  roots  only  in  the  open  left  half 
plane.  We  could  use  the  Hurwitz  criterion  (cf.  Appendix  A)  to  obtain  a 
system  of  nonlinear  inequalities  which  describes  the  feasible  region.  A 
simpler  way  is  described  as  follows:  Consider  the  region 

X  =  (xeKk-1  |  g(x,w,j)  <  0  for  all  we[0,o°),  j=l,2}  (5.3) 

where 

g(x,w,l )  =  -  |r(x,iw)  |2 
g(x,w,2)  =  -  |s(x,iw)  |2 

It  is  obvious  that  both  g(x,w,j)  and  3g(x,w,j)/8x  are  continuous  in 
Kk-1x[0,°°)x{l  ,2).  The  set  X  contains  all  x  such  that  r(x,z)  and  s(x,z) 
have  no  roots  on  the  imaginary  axis  (in  the  complex  z-space).  Thus  X  con¬ 
sists  of  a  number  of  disconnected  open  sets  each  of  which  is  either  entire¬ 
ly  feasible  or  not  feasible.  If  the  initial  trial  x°  lies  in  a  feasible 
component,  so  will  all  subsequent  trials  if  the  step  selection  rule  does 
not  take  too  large  a  step  landing  into  an  infeasible  component.  To  prevent 
that,  we  only  accept  those  points  x  such  that  r(x,z)  and  s(x,z)  have  roots 
all  in  the  open  left  half  z-plane  and  reject  all  other  points.  In  other 
words,  the  feasible  region  is  implicitly  described  by  X  and  the  above  test, 
but  only  the  function  g(x,w,j)  is  needed  in  the  determination  of  feasible 
descent  directions. 

The  objective  function  for  the  minimax  formulation  is 

Re  r(x,iw)s(x,-iw) 

f(x,w)  =  - 

|Im  r(x,iw)s(x,-iw)  | 

whose  value  at  (x,w)  is  equal  to  -  cot  |  Arg  r(x,iw)/s(x,iw) | .  Although  the 
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negative  part  of  f(x,w)  is  unbounded,  it  is  never  needed  in  determining  the 
maximum  of  f(x,w)  over  we[0,°°).  From  the  fact  that  as  w  -*■  °°, 

f(x,w)  ->0 

g(x,w,j )  +  -  «  j  =  1,2 

it  is  not  difficult  to  show  that  the  results  of  Chapter  IV  still  apply  so 
long  as  the  positive  part  of  f(x,w)  is  bounded,  which  is  true  for  those  x 
representing  methods  which  are  A(Q)-stable.  Hence,  if  we  can  find  an  x^eX 
which  represents  an  A(0)-stable  A^-stable  method  satisfying  the  strict  root 
condition,  then  we  can  use  the  M-algorithm  to  find  the  infimum  of  the  func¬ 
tion 

<j)(x)  =  max  f(x,w) 
we[0,°°) 


over  the  set 


X{x°)  =  (xeX  |  (x )  £  <J> (x°)  > 

together  with  the  test  on  r(x,z)  and  s(x,z)  mentioned  earlier.  The  M- 
algorithm  will  then  either  conclude  that  no  stationary  point  exists  or  give 
a  point  x*  such  that  the  value  of  <j>(x*)  is  close  to  the  value  of  a  local 
infimum  of  <j>  on  X(x°).  Note  that  since  X(x°)  is  not  closed,  the  infimum 
may  not  be  attained  by  any  point  in  X(x°). 

The  initial  trial  x°  can  be  found  by  first  finding  an  A(0)-stable 
method  and  then  by  varying  some  of  the  parameters  bp,b1,...,bk_1  such  that 
the  magnitude  of  the  error  constant  is  without  violating  the  feasibility 
of  the  new  point  in  the  x-space. 


5.3  NUMERICAL  RESULTS 

The  M-algorithm  described  in  section  4.3  is  implemented  to  run  on 
the  Cyber  175  using  double  precision  arithmetic.  The  numerical  results 
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will  be  given  after  a  brief  discussion  on  the  Fortran  program. 

The  evaluation  of  <f>(x)  is  in  two  steps.  First,  we  test  if  the 
method  represented  by  x  is  A(0)-stable.  From  (3.2),  the  method,  being 
Actable  (see  constraints  given  in  the  previous  section),  is  A(0)-stable 
if  and  only  if  the  following  condition  is  not  satisfied  by  any  positive 
real  w  of  Im  r(x,iw)s(x,-iw)  =  0: 

Re  r(x,iw)s(x,-iw)  <  0 

The  test  involves  finding  the  positive  roots  of  an  even  polynomial  of  de¬ 
gree  2 ( k-1 )  since 


r(x,iw)s(x,-iw)  = 


£  (-l)jC(2j)w2j  +  iw  V  (-1 )  J+1C(2j+l  )w2j 


u(x,w)  +  iv(x,w) 


where 


C(j)  =  e  (-1  )va  b 


V  j-v 


for  j  =  0,1 ,..,,2k-l 


Note  that  from  (1.9),  it  can  be  proved  that  C(j)  =  0  for  j  >_ k,  j  even.  If 

the  method  represented  by  x  is  A(0)-stable,  then  <f>(x)  is  found  by  locating 

all  the  local  maxima  of  f(x,w)  on  0  <  w  <  ®.  Again,  we  have  to  solve  for 
the  positive  roots  of  an  even  polynomial  of  degree  at  most  3( k-1 )  since  it 

can  be  shown  that  for  w  >  0  such  that  Im  r(x,iw)s(x,-iw)  f  0, 

3f  v(x,w)[3u(xpw)/3w]  -  u(x,w)[3v(x,w)/3w] 


3f 

- (x,w)  =  - 

3w 


|v(x,w)  |v(x,w) 


where  the  numerator  simplifies  to 


l_3(k-l  )/2j  v 

£  H)V  z  (2v-4j+l )C (2 j )C(2v-2j+l )w^ 
v=0  j=0 


The  set  of  local  maxima  contains  R(x)  as  a  subset. 
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Similarly,  the  set  Q(x)  is  contained  in  the  set  of  local  maxima 

over  we[Q,«)  of  the  functions  g(x,w,j),  j  =  1,2,  given  by  (cf.  (5.3)) 

k-1  .•  pi 

g(x,w,l )  =  -  s  (-1  rA(2j)v/  J 
j=0 

g(x,w,2)  =  -  Z  (-l)JB(2j)w2j 
0=0 

where  for  j  =  0,1 . 2k, 

a(ji=  i 

v=Q 

B(J)  =  i  H)\bj.v 

v=0  J 

As  remarked  in  section  4.3,  the  sets  R£(x)  and  Qy(x)  are  discre¬ 
tized  so  that  ||z  ( x ) | I  can  be  approximated  by  the  solution  of  a  convex 

quadratic  programming  problem.  The  discretization  contains,  in  addition  to 
points  in  R(x)  and  Q(x),  points  of  the  form  w±<5  in  R£(x)  for  each  weR(x) 
and  (w±S,j.)  in  Q^x)  for  each  (w,j)eQ(x).  Such  points  are  obtained  using 
the  following  procedure:  Let  weR(x).  Start  with  an  initial  value  for  6 
given  by  *^e".  If  f(x,w)  -  f(x,w -5)  e  then  w±<5  are  the  points  chosen. 

Otherwise,  set  6  =  <$£/[f(x,w)  -  f(x,w-<5)]  and  repeat  the  above  test.  The 
procedure  for  obtaining  (wt6,j)  in  Qy(x)  is  similar:  Let  (w,j)eQ(x).  If 
g(x,w,j)  -  -  y,  set  6  -  0.  Suppose  that  g(x,w,j)  >  -  y.  Start  with  6  = 

47.  If  g(x,w-<$,j)  >_  -  y  then  (w±6,j)  are  accepted;  else  set  6  =  6[g(x,w,j) 
+  y]/[g(x,w,j)  -  g(x,w-6 ,j )]  and  repeat.  Empirical  results  show  that  the 
convergence  rate  is  faster  if  we  include  the  additional  test  to  check  if 
(0,j)  is  in  Qp(x). 

For  given  x  and  w,  the  vector  9f(x,w)/3x  can  be  found  using  two 
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linear  recurrence  relations,  the  derivation  of  which  is  straightforward: 


3f 


3X 


(x.w)  =  - 


1 


T,R  +  T,I 

1  v  2  v 


V  —1,3,... 


1 


I  In*  r(x,iw)s(x,-iw)  |  |  -  T3(wlv_1+— )  v  =2,4, 


where 

x  ,  1  <  v  <  k-1 ,  is  the  v-th  component  of  x, 

T-j  =  2 [Re  s(x,iw)  +  F(Im  s(x,iw))] 

T2  =  Im  r(x,iw)  +  F(Re  r{x,iw)) 

T3  =  Re  r(x,iw)  -  F(Im  r{x,iw)) 

T4  =  2w[Im  s(x,iw)  -  F(Re  s(x,iw))] 

u(x,w) 

F  =  - 

v(x,w) 

and  Rv,  Iv  are  defined  by  the  recurrence  relations 

1  2 

R-i-l  Rv  -  ■  w  R^  ,  v  =  3,5,... 

v  ~c 

1 1  -  w  I\>  ~  ”  w  v  =  ^»^*‘** 


Similarly,  the  vector  3g(x,w,j)/3x  for  given  (x,w,j)  can  be  ob¬ 
tained  using  the  above  recurrence  relations: 


4  Re  r(x,iw)  Ry 

j  =  1 

v  = 

1,3,... 

9g 

—  (x,w,j)  =  ( 

-  4  w  Im  r(x,iw)  R  , 
v- 1 

j  =  1 

V  = 

2,4,... 

3X 

V 

-  2  Im  s(x,iw)  I 

j  =  2 

V  = 

1,3,... 

^2  Re  s(x,iw)  [wly_1  +  - ] 

v+1 

CM 

II 

V  = 

2,4,... 

The  program  starts  with  e-j  =  y-j  =  =  .001  and  decreases  the 

parameters  by  a  factor  of  1/2  in  each  iteration,  i.e., 

ek+l  =  ek/2  pk+l  =  pk/2  pk+l  =  pk/2 
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In  case  the  points  generated  by  the  (e,y,p)-algorithm  with  e  =  ek,  p  =  vk. 
p  =  Pk  converge  to  a  nonstationary  point,  the  program  makes  one  more  at¬ 
tempt  by  setting  ek+1  =  Pk+1  =  Pk+1  =  .1  before  an  abnormal  end  is  signal¬ 
led. 

Before  we  can  use  the  M-algorithm  to  find  A*(a),  an  initial 
feasible  point  b°  =  b°[A]  has  to  be  known.  It  can  be  supplied  by  the  user 
as  an  input  to  the  program.  The  program  also  has  the  capability  of  finding 
b°[A]  if  it  is  given  a  point  b°[A']  representing  a  method  with  error  con¬ 
stant  of  magnitude  ( A ' ) k  instead  with  A'  4  A.  The  program  first  applies 
the  M-algorithm  for  the  value  A'  to  obtain  a  point  b*  =  b*[A'].  Then  it 
tests  if  the  vector  b'  whose  v-th  component  is  given  by 


K 

v  =  1,3 

A  is 

(-)kb* 

v  -  0,2 

A 


is  feasible  for  the  value  A  and  represents  an  A(0)-stable  method.  Note 

that  from  the  definition  of  b’,  the  corresponding  method  has  error  constant 

of  magnitude  (A)k.  The  above  definition  for  b'  works  satisfactorily  for  odd 

k  and  for  even  k  when  A  >  A'  in  the  sense  that  in  most  cases,  b'  is  a 

feasible  point  and  does  represent  an  A(0)-stable  method.  For  the  case  when 

k  is  even  and  A'  >  A,  we  make  use  of  an  empirical  observation  that  b* 

always  tends  to  zero  and  hence  the  polynomial  s(b*,z)  must  have  at  least 

one  real  (negative)  root.  Instead  of  fixing  A  and  trying  to  find  b'  so 

k 

that  the  corresponding  method  has  error  constant  of  magnitude  (a)  ,  we  de¬ 
crease  the  largest  negative  real  root  of  s(b*,z),  say  y,  by  a  factor  of  .1 
and  define  the  new  polynomial  s(b',z)  by 
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z  -  .ly 

s(b',z)  =  - s(b*,z) 


z  -  y 

It  is  straightforward  to  show  that  for  v  =  l,2,...,k-l, 


bv  = 


bv* 


+  .  9y 


k 

E 

j=v+l 


with  bg  =  b^  =  0  and  that  the  value  A  corresponding  to  b'  is  smaller  than 
A1.  The  above  definitions  for  b1  are  only  heuristics  for  finding  an  initi¬ 
al  feasible  A(0)-stable  method  for  the  M-algorithm.  In  all  cases,  if  the 
vector  b'  is  not  feasible  or  does  not  represent  an  A(0)-stable  method,  the 
program  will  try  to  proceed  in  smaller  steps. 

The  positive  real  roots  of  a  given  polynomial  are  found  by  first 
generating  the  Sturm  sequence  for  the  polynomial  and  its  derivative  in  the 
interval  [O,00)  and  then  using  the  ZEROIN  routine  (Shampine  and  Allen,  1973, 


pp. 244-246)  which  locates  a  root  of  a  continuous  function  inside  a  given 
interval  by  a  combination  of  bisection  and  secant  rules.  The  convex  qua¬ 
dratic  programming  problem  is  solved  using  the  first  symmetric  variant  of 
the  simplex  method  (Van  De  Panne,  1975,  pp. 270-280). 

Numerical  results  for  A*(a)  for  a  finite  number  of  points  on  0  < 
A  <  00  for  the  cases  when  k  =  4, 5, 6, 7  are  plotted  in  figures  5.5,  5.6,  5.7, 
5.8,  respectively.  The  isolated  symbols  +  represent  the  (A,a)  values  with 
A  =  |ck+i  |  and  a  being  the  angle  of  absolute  stability  for  the  following 
methods:  BDF  (Backward  Differentiation  Formula),  CHEB2,  CHEB3,  CHEB4  (Gup¬ 
ta,  1976),  and  I£,  (Gupta  and  Wallace,  1975).  The  BDF  of  order  6  and 
the  other  methods  in  Gupta  (1976)  lie  outside  the  region  described  in  the 
plots  and  hence  are  not  included.  The  table  following  each  figure  contains 
the  values  of  the  parameters  bp,^,. ..  ,bk_i  for  the  method  corresponding  to 
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b0 

bl 

b2 

b3 

.0 

.0022 

.4165 

.6103 

.0 

.0041 

.4526 

.6241 

.0 

.0066 

.4913 

.6389 

.0 

.0094 

.5284 

.6530 

.0 

.0136 

.5787 

.6716 

.0 

.0198 

.6411 

.6953 

.0 

.0282 

.7159 

.7226 

.0 

.0383 

.7963 

.7513 

.0 

.0462 

.8544 

.7714 

.0 

.0597 

.9458 

.8025 

.0 

.0950 

1.1544 

.8700 

.0 

.1164 

1.2675 

.9048 

.0 

.1428 

1.3980 

.9434 

.0 

.2021 

1 . 6630 

1.0177 

.0 

.2864 

1.9979 

1.1049 

.0 

.4367 

2.5206 

1.2290 

.0 

.6673 

3.2156 

1.3773 

.0 

1.2635 

4.7096 

1.6522 

.0 

3.5655 

9.0 

2.2637 

.0 

13.2348 

21.0 

3.4392 

Table  5.1  k  =  4 
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b0 

bl 

b2 

b3 

b4 

.0 

1.366 

.816 

3.238 

.373 

.0 

1.366 

.817 

3.238 

.373 

.0 

1.366 

.819 

3.238 

.374 

.0 

1.366 

.834 

3.238 

.381 

.0 

1.366 

.863 

3.238 

.394 

.0 

1.366 

.981 

3.238 

.448 

.0 

1.675 

1.998 

3.568 

.837 

.0 

3.081 

3.359 

4.709 

1.068 

.0 

6.553 

6.171 

6.733 

1.381 

.0 

10.660 

9.036 

8.518 

1.606 

.0 

15.241 

11.929 

10.141 

1.786 

.0 

20.226 

14.837 

11.648 

1.938 

.0 

31.206 

20.687 

14.417 

2.189 

.0 

37.131 

23.622 

15.707 

2.296 

.0 

49.725 

29.508 

18.144 

2.487 

.0 

70.236 

38.363 

21.527 

2.729 

.0 

82.374 

43.292 

23.296 

2.847 

.0 

113.365 

55.141 

27.296 

3.099 

.0 

146.699 

67.011 

31.025 

3.316 

.0 

182.088 

78.895 

34.544 

3.508 

Table 

5.2  k  = 

5 
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b0 

bl 

b2 

b3 

b4 

b5 

.0 

.233 

3.470 

3.579 

4.218 

1.295 

.0 

.298 

3.928 

3.882 

4.459 

1.329 

.0 

.385 

4.491 

4.244 

4.738 

1.367 

.0 

.500 

5.188 

4.678 

5.061 

1.410 

.0 

.653 

6.054 

5.200 

5.436 

1.459 

.0 

.749 

6.565 

5.500 

.5.645 

1.486 

.0 

.860 

7.137 

5.829 

5.869 

1.514 

.0 

1.139 

8.496 

6.586 

6.371 

1.576 

.0 

1.314 

9.302 

7.021 

6.650 

1.609 

.0 

1.756 

11.227 

8.026 

7.273 

1.682 

.0 

2.034 

12.375 

8.605 

7.620 

1.721 

.0 

2.741 

15.130 

9.943 

8.392 

1.806 

.0 

3.188 

16.779 

10.716 

8.821 

1.852 

.0 

4.713 

22.032 

13.061 

10.065 

1.978 

.0 

6.082 

26.405 

14.907 

10.992 

2.068 

.0 

7.871 

31.780 

17.076 

12.033 

2.164 

.0 

12.836 

45.421 

22.209 

14.334 

2.363 

.0 

17.215 

56.451 

26.077 

15.951 

2.494 

.0 

26.848 

78.745 

33.367 

18.794 

2.709 

.0 

37.431 

101.254 

40.213 

21.280 

2.884 

Table  5.3  k  =  6 
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b0 

V 

b2 

b3 

b4 

b5 

b6 

.0 

9.750 

26.031 

29.627 

15.260 

10.578 

1.898 

.0 

14.853 

35.990 

37.941 

18.553 

11.942 

2.048 

.0 

46.134 

87.855 

75.803 

31.785 

16.798 

2.506 

.0 

133.735 

206.759 

148.501 

53.244 

23.442 

3.020 

.0 

151.480 

228.727 

160.819 

56.584 

24.388 

3.086 

.0 

178.912 

261.810 

178.928 

61.381 

25.716 

3.176 

.0 

202.449 

289.483 

193.719 

65.211 

26.751 

3.245 

.0 

226.560 

317.238 

208.269 

68.907 

27.732 

3.308 

.0 

251.198 

345.065 

222.599 

72.485 

28.665 

3.368 

.0 

276.330 

372.957 

236.729 

75.959 

29.556 

3.424 

.0 

301.921 

400.908 

250.674 

79.337 

30.410 

3.477 

.0 

327.951 

428.912 

264.450 

82.628 

31.230 

3.526 

.0 

354.396 

456.964 

278.068 

85.840 

32.020 

3.574 

.0 

381.237 

485.061 

291.539 

88.980 

32.783 

3.619 

.0 

436.028 

541.376 

318.076 

95.061 

34.236 

3.703 

Table  5.4  k  =  7 
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the  symbol  H  in  ascending  order  of  A.  The  numerical  values  for  (a,A*(a)) 
and  the  corresponding  (rG)1/,k  and  f  are  listed  in  Appendix  B. 

5.4  CONCLUSION 

From  the  empirical  results,  it  is  observed  that  the  parameter  b^ 
corresponding  to  the  method  found  by  the  M-algorithm  is  always  zero  (or 
very  close  to  zero)  which  implies  that  a (c)  has  a  root  at  ?  =  -1.  Moreo¬ 
ver,  except  for  the  root  s  =  1  of  p(?)  and  the  root  c  =  -1  of  a(c),  the 
other  roots  of  p(c)  and  ct(?)  are  close  together  and  they  seem  to  be  ap¬ 
proaching  the  point  ?  *  1  as  A  •»■  »,  It  is  interesting  to  note  that  for  all 
the  methods  found,  the  ratio  (rG)1/k  to  |ck+1|1/k  is  not  very  large  (<  2), 
and  that  f  <  1.3.  But  each  of  the  methods  has  a  pair  of  roots  of  p(e)  of 
magnitude  close  to  1.  Thus,  the  global  error  accumulated  when  any  of  the 
methods  is  used  is  not  that  much  affected  by  the  positions  of  the  roots  of 
p(0  provided  a  fixed  stepsize  is  used  throughout. 

We  could  compare  the  BDF  with  the  methods  found  by  the  M- 
algorithm  as  follows:  Since  the  k-th  order  k-step  BDF  has  error  constant  of 
magnitude  l/(k+l),  we  consider  the  ratio 

\  '  aBDF[k] 
f  -  A*((l/(k+l))1/k) 

where  aBDp[|<]  is  the  angle  of  absolute  stability  associated  with  the  BDF. 
For  4  £  k  i  6  (the  7-th  order  7-step  BDF  is  not  A(O)-stable) ,  the  above  ra¬ 
tio  is  greater  than  3.  In  other  words,  the  method  found  by  the  M-algorithm 
having  error  constant  of  magnitude  l/(k+l)  is  at  least  3  times  "closer  to 
being  A-stable"  than  the  BDF.  On  the  other  hand,  for  4  <  k  <.  6,  the  ratio 

(l/(k+l))1/k 
Ak 


>  1.3 
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where  Ak  is  the  value  of  a  such  that  A*(a^)  =  aBDF^kj.  Thus  for  the  method 
corresponding  to  ( a k , A* ( a k ) )  which  has  the  same  angle  of  absolute  stability 
as  the  BDF,  the  error  constant  has  magnitude  less  than  (10/13) k  times  that 
of  the  BDF.  This  implies  that  if,  in  an  ODE  solver,  the  steps ize  is  chosen 
to  be  proportional  to  the  k-th  root  of  the  reciprocal  of  the  magnitude  of 
the  error  constant,  then  the  method  corresponding  to  (Ak,A*(Ak))  is  1.3 
times  more  efficient  than  the  BDF.  However,  numerical  tests  (based  on 
DIFSUB  (Gear,  1971)  with  appropriate  modifications)  show  that  for  problems 
which  can  be  successfully  handled  by  the  BDF,  the  method  corresponding  to 
(Ak,A*(Ak))  is  only  about  .7  times  as  efficient  as  the  BDF.  The  reason  may 
be  that  our  theory  did  not  take  into  account  the  local  instabilities  that 
could  arise  from  step  changing  and  order  changing.  This  is  especially 
likely  since  the  p(^)  and  a(c)  polynomials  associated  with  the  methods  have 
roots  which  are  much  nearer  the  unit  circle  than  those  of  the  BDF. 

A  final  point  is  that  we  can  see  from  the  plots  that  some  of  the 
methods  obtained  by  Gupta  and  Wallace  (1975)  and  Gupta  (1976)  are  close  to 
the  (a,A*(a))  values. 
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POLYNOMIALS  HAVING  ROOTS  IN  THE  OPEN  LEFT  HALF  COMPLEX  PLANE 
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The  following  is  a  list  of  equivalent  conditions  each  of  which  is 
necessary  and  sufficient  for  a  polynomial 

p(z)  =  cnzn  +  cn_-|zn  1  +  ...  +  C-|Z  +  c0  cn  >  0 

to  be  Hurwitzian,  i.e.,  to  have  roots  all  in  the  open  left  half  complex 
plane: 

(1)  (Hurwitz  criterion)  The  leading  principal  minors  of  the  nxn 
matrix  H  whose  (i,j)-th  element  is  given  by 

Hij  =  cn-i-2j  1*J  = 

are  positive  (Krall,  1967,  p.48).  A  stronger  statement  (Lienard  and  Chi- 
part  criterion)  is  in  Gantmacher  (1959,  p. 221 ) . 

(2)  (Routh  criterion)  The  elements  in  the  first  column  of  the 
tableau  R  are  positive,  where  the  elements  in  R  are  generated  by  the  fol¬ 
lowing  scheme  (Krall,  1967,  p. 52) :  Start  with  the  first  two  rows  whose 
elements  are 

cn  cn-2  cn-4  cn-6  ••• 

cn-l  cn-3  cn-5  cn-7  *** 

Every  row  after  the  first  and  second  in  the  tableau  is  determined  by  the 
two  preceding  rows  according  to  the  following  rule:  If  the  elements  of  the 
two  preceding  rows  are 

9o  9]  92  93 

h0  hl  h2  h3 

then  the  elements  q0,q-|,...  of  the  current  row  are  given  by 

g0 . 

qj "  gj+i  "  r^j+i 
nQ 


j  =  0,1,... 
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The  tableau  consists  of  the  first  n+1  rows  generated  by  the  scheme. 

(3)  (L.  Cremer- Leon hard  separation  criterion)  The  zeros  {wy}  and 

(wv)  of  the  polynomials  Re  p(if'-w)  and  Im  pCif-w)//^,  respectively,  sa¬ 
tisfy 

0  >  w|  >  w-]  >  w£  >  w2  >  ... 


(Lehnigk,  1966,  p. 1 54) . 

(4)  The  polynomial 

2  n-1  ,  ■  .  ,  n-2  n-3 

cn-l2  +  ^cn-lcn-2“cncn-3^z  +  cn-lcn-3z  + 

+  (cn-lcn-4-cncn-5)zn"4  +  '•••  + 


+  r  r  ,n-2m+l  ,  /  .  n-2m 

cn-lcn-2m+lz  +  (cn-lcn-2rrfcncn-2m-l)z  + 


is  Hurwitzian  (Krall,  1967,  p.47). 

We  next  give  some  necessary  conditions  for  p(z)  to  be  Hurwitzian. 
Suppose  that  p(z)  is  Hurwitzian.  Then  it  is  obvious  that 

Cj  >  0  j  =  0,1,..., n 

Using  (4),  it  follows  that 

cn-lcn-j  "  cncn-j-l  >  0  ^  2»4»“* 

Moreover,  from  (3),  using  the  Dougall's  inequality  for  symmetric  functions 
(Mitrinovic,  1970,  p.97),  it  can  be  shown  that  the  following  inequalities 
are  satisfied: 


2 

C2 

4 

-  s 

c4 

n-2 

s.  -  -  -  •*. 

cn-2  >  n 

cn 

n 

c0 

—  n-2 

c2 

-  4 

cn-4  "  2 

cn-2 

2 

C3 

4 

C5 

n-4 

cn-3 

2  cn 
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when  n  is  even,  and 


2 

c2 

> 

4 

C4 

—  >  .  • . 

n-3 

>  - 

cn-3 

n-1 

> - 

cn-l 

n-1 

co 

n-3 

C2~ 

“  4 

cn-5 

“  2 

c  n-3 

2 

C3 

> 

4 

C5 

—  >  •  • « 

n-3 

cn-2 

n-1 

>  - 

cn 

n-1 

C1 

n-3 

C3~ 

-T- 

cn-4 

“  2 

cn-2 

when  n  is  odd 
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APPENDIX  B 


LIST  OF  A,  A*(A ),  (rG)1/k,  f  VALUES  FOR  METHODS 


FOUND  USING  THE  M-ALGORITHM 
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A 

A*  (A) 

(rG)1/k 

f 

.3815 

12.17 

.5069 

1.1388 

.3848 

16.45 

.5086 

1.1381 

.3883 

20.55 

.5103 

1.1374 

.3916 

24.06 

.5118 

1.1366 

.3959 

28.24 

.5137 

1.1353 

.4010 

32.87 

.5158 

1.1336 

.4069 

37.53 

.5182 

1.1314 

.4130 

41.78 

.5205 

1.1288 

.4172 

44.43 

.5221 

1.1268 

.4236 

48.08 

.5243 

1.1236 

.4372 

54.59 

.5287 

1.1162 

.4441 

57.34 

.5307 

1.1123 

.4517 

60.03 

.5328 

1.1078 

.4660 

64.32 

.5364 

1.0993 

.4823 

68.27 

.5450 

1.0896 

.5050 

72.47 

.5692 

1.0769 

.5310 

76.06 

.5919 

1.0636 

.5767 

80.32 

.6403 

1.0445 

.6687 

84.85 

.7323 

1.0204 

.8190 

87.77 

.8888 

1.0098 

Table  B.l  k  =  4 


A 

A*  (a) 

,4045 

.01 

,4046 

.07 

,4048 

.20 

,4063 

1.22 

,4091 

3.21 

,4197 

10.47 

,4821 

43.49 

5296 

57.68  ’ 

5923 

68.83 

6361 

73.72 

6704 

76.54 

6988 

78.40 

7448 

80.75 

7641 

81.54 

7977 

82.73 

8394 

83.91 

8594 

84.38 

9010 

85.23 

9361 

85.81 

9666 

86.25 

Table  B.2 

(rG) 

f 

.6233 

1.2074 

.6233 

1.2073 

.6232 

1.2072 

.6228 

1.2066 

.6218 

1.2052 

.6184 

1.2001 

.6238 

1.1649 

.6602 

1.1321 

.7300 

1.0937 

.7776 

1.0716 

.8206 

1.0606 

.8517 

1.0563 

.9145 

1.0484 

.9385 

1.0451 

.9847 

1.0393 

1.0417 

1.0326 

1.0708 

1.0296 

1.1315 

1.0256 

1.1848 

1.0231 

1.2300 

1.0210 

k  =  5 

A 

A*  (a) 

5677 

40.52 

5763 

43.54 

5858 

46.59 

5965 

49.67 

6085 

52.72 

6149 

54.23 

6217 

55.72 

6363 

58.64 

6442 

60.06 

6610 

62.80 

,6700 

64.12 

,6893 

66.65 

.6995 

67.85 

.7230 

70.29 

.7473 

72.42 

.7681 

74.00 

.8106 

76.66 

.8379 

78.06 

.8821 

79.91 

.9174 

81.12 

Table  B.3 

(rg)l/k 

f 

.7072 

1.1858 

.7139 

1.1807 

.7234 

1.1750 

.7280 

1.1687 

.7412 

1.1617 

.7467 

1.1579 

.7548 

1.1540 

.7662 

1.1457 

.7694 

1.1414 

.7915 

1.1324 

.7986 

1.1278 

.8169 

1.1182 

.8291 

1.1134 

.8492 

1.1031 

.8740 

1.0930 

.8985 

1.0852 

.9420 

1.0711 

.9716 

1.0632 

1.0166 

1.0522 

1.0570 

1.0448 

k  =  6 

A 

A*  (A) 

(TG)1/k 

A 

r 

.7131 

54.28 

.9039 

1.1676 

.7430 

58.98 

.9322 

1.1524 

.8343 

68.85 

1.0306 

1.1122 

.9351 

75.26 

1.1381 

1.0785 

.9479 

75.87 

1.1545 

1.0749 

.9653 

76.64 

1.1722 

1.0703 

.9786 

77.19 

1.1875 

1.0669 

.9908 

77.67 

1.2021 

1.0640 

1.0022 

78.09 

1.2134 

1.0613 

1.0129 

78.47 

1.2270 

1.0589 

1.0229 

78.81 

1.2374 

1.0568 

1.0324 

79.12 

1.2483 

1.0557 

1.0414 

79.41 

1.2595 

1.0548 

1.0499 

79.66 

1.2674 

1.0540 

1.0658 

80.12 

1.2865 

1.0524 

Table  B.4 

k  =  7 
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